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Under the driving conditions, a typical geared system may be subjected to large 
dynamic loads. Also, the vibration level of the geared system is directly related to the 
noise radiated from the gear box. Accordingly, a good understanding of the steady 
state dynamic behavior of the system is required in order to design reliable and quiet 
transmissions. It is the main focus of this study with the scope limited to a system 
containing a spur gear pair with backlash and periodically time-varying mesh stiffness, 
and rolling element bearings with clearance type non-linearities. The internal static 
transmission error excitation at the gear mesh, which is of importance from high 
frequency noise and vibration control view point, is considered in the formulation in 
sinusoidal or periodic form. 

A dynamic finite element model of the linear time-invariant (LTI) system is 
developed. Effects of several system parameters, such as torsional and transverse 
flexibilities of the shafts and prime mover/load inertias, on free and forced vibration 
characteristics are investigated. Several reduced order LTI models are developed and 
validated by comparing their eigen solutions with the finite element model results. 
Using the reduced order formulations, a three-degree of freedom dynamic model is 
developed which includes non-linearities associated with radial clearances in the radial 
rolling element bearings, backlash between a spur gear pair and periodically varying 
gear meshing stiffness. As a limiting case, a single degree of freedom model of the 
spur gear pair with backlash is considered and mathematical conditions for tooth 
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separation and back collision are defined. Both digital simulation technique and 
analytical methods such as method of harmonic balance and the method of multiple 
scales have been used to develop the steady state frequency response characteristics for 
various non-linear and/or time-varying cases. Difficulties associated with the 
determination of the multiple solutions at a given frequency in the digital simulation 
technique have been resolved. The proposed formulation has been validated by 
comparing the predictions with the results of two benchmark experiments reported in 
the literature. Several key system parameters such as mean load and damping ratio are 
identified and their effects on the non-linear frequency response are evaluated 
quantitatively. Other fundamental issues such as the dynamic coupling between non- 
linear modes, dynamic interactions between component non-linearities and time-varying 
mesh stiffness, and the existence of subharmonic and chaotic solutions including routes 
to chaos have also been examined in depth. 
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CHAPTER I 


INTRODUCTION 
1.1. PROBLEM FORMULATION 

Dynamic analysis of geared systems is an essential step in design due to two 
reasons. First, under the driving conditions, a typical geared system is subject to 
dynamic forces which can be large. Therefore, the prediction of dynamic loads, 
motions or stresses is needed in developing reliable gear trains. Second, the vibration 
level of the geared system is directly related to the noise radiated from the gear box. An 
attempt in designing quiet gears requires a good understanding of the dynamic behavior 
of the system and the gear mesh source. Accordingly, the main objective of this study 
is to develop accurate mathematical models of a generic geared rotor-bearing system 
shown in Figure 1.1a. Of interest here is to investigate several key modelling issues 
which have not been addressed in the literature, such as system non-linearities and 
time-varying mesh stiffness. 

The generic geared system shown in Figure 1.1a consists of a single spur gear 
mesh of ratio v g =dg 2 /dgi > rolling element bearings, a prime mover driving the system 

at Q s j speed and a typical inertial load. The system also includes other elements such 

as couplings and flywheel. A discrete model of the system is shown in Figure 1.1b. 
Here, shafts are represented by discrete translational springs k sl and k s2 , translational 
dampers c sl and c s2 , and torsional springs K sl and K $2 . The gear mesh is represented 
by a time-varying mesh stiffness k^(t) and a non-linear displacement function f h 
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Ig2i m g2,dg2 

(b) 


Figure 1.1. a) A generic geared rotating system, b) discrete model of geared rotating 



which includes gear backlash. Further, linear time-invariant (LTI) mesh damping c^ is 

considered here. The rolling element bearings are defined by a time-invariant radial 
stiffness k^ subject to a non-linear displacement function fj,, and an LTI damping 
coefficient c^. The prime mover and load are modeled as purely torsional elements of 
inertias Ip and 1^, respectively. The mean rotational speeds and ft s 2 and the 

geometric end conditions are such that gyroscopic effects are not seen. 

The gener aliz ed displacement vector (q(t)}, associated with the inertia elements, 

consists of angular displacements 8 and transverse displacements x and y . The 
governing equation of motion for the non-linear, time-varying multi-degree of freedom 
model can be given in the general form as 

[M]{q"(t)} + [C]{q'(t )} + [K(t)]{f (q(t))} = {F(t)} (1.1) 


where [m] is the time-invariant mass matrix and {q(t)} * s the displacement vector. 
Here, damping matrix [c] is assumed to be LTI type, as the effect of the tooth 

separation and time-varying mesh properties on mesh damping are considered 

negligible; validity of this assumption will be examined later. The stiffness matrix 
[K(t>] is considered to be time-varying, given by a periodically time varying matrix 

[K(t)j = [K(t + 2 ti / Q h )] where is the fundamental gear mesh frequency. The 

non-linear displacement vector {f (q(t ))} includes the radial clearances in bearings and 
the gear backlash, and the forcing vector |F(t)J consists of both external excitations 

due to torque fluctuations, mass unbalances and geometric eccentricities, and an 
internal static transmission error excitation. 
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1.2. OBJECTIVES 

Specific objectives of this study are given as follows; each chapter, written in the 
journal paper style, deals with one major objective. 

First, a dynamic finite element model of the linear time-invariant (LTT) system 
given in Figure 1.1a is developed. Effects of several system parameters, such as 
torsional and transverse flexibilities of the shafts and prime mover/load inertias, on free 
and forced vibration characteristics are investigated. Three different reduced order LTI 
models will be derived and the conditions under which these are valid will be 
determined by comparing the eigen solutions with the finite element model results. 
Development and verification of such a reduced order (with a very few degrees of 
freedom) linear model is an essential step before the non-linear dynamic behavior is 
analyzed [Chapter I]. 

Second, non-linear frequency response characteristics of a spur gear pair with 
backlash are examined for both external and internal excitations. The internal excitation 
is of importance from the high frequency noise and vibration control viewpoint and it 
represents the overall kinematic or static transmission error. Such problems may be 
significantly different from the rattle problems associated with external, low frequency 
torque excitation. Two solution methods, namely the digital simulation technique and 
the method of harmonic balance have been used to develop the steady state solutions 
for the internal sinusoidal excitation. Difficulties associated with the determination of 
the multiple solutions at a given frequency in the digital simulation technique have been 
resolved as one must search the entire initial conditions map. Such solutions and the 
transition frequencies for various impact situations are found by the method of 
harmonic balance. Further, the principle of superposition can be employed to analyze 
the periodic transmission error excitation and/or combined excitation problems 
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provided the excitation frequencies are sufficiendy apart from each other. Predictions 
are compared with the limited experimental data available in the literature [Chapter II]. 

Third, non-linear frequency response characteristics of a geared rotor-bearing 
system are examined. A three degree of freedom dynamic model is developed which 
includes non-linearities associated with radial clearances in the radial rolling element 
bearings and backlash between a spur gear pair; time-invariant gear meshing stiffness is 
assumed. The bearing non-linear stiffness function is approximated for convenience 
sake by a simple model which is identical to that used for the gear mesh. This 
approximate bearing model has been verified by comparing steady state frequency 
spectra. The applicability of both analytical and numerical solution techniques to the 
multi degree of freedom non-linear problem is investigated. Proposed theory is 
validated by comparing the results with available experimental data. Several key 
issues such as non-linear modal interactions and differences between internal static 
transmission error excitation and external torque excitation are discussed. Additionally, 
parametric studies are performed to understand the effect of system parameters such as 
bearing stiffness to gear mesh stiffness ratio, alternating to mean force ratio and radial 
bearing preload to mean force ratio on the non-linear dynamic behavior. A criterion 
used to classify the steady state solutions is presented and the conditions for chaotic, 
quasi-periodic and subharmonic steady state solutions are determined. Two typical 
routes to chaos observed in this geared system are also identified [Chapter III]. 

Fourth, this study extends the non-linear single degree of freedom spur gear pair 
model of Chapter II and multi-degree of freedom geared rotor-bearing system model of 
Chapter IE by including time-varying mesh stiffness k|,(t), and investigates the effect 

of kh(t) on the frequency response of lightly and heavily loaded geared systems. 
Interactions between mesh stiffness variation and system non-linearities associated with 
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gear backlash and radial clearances in rolling element bearings are also considered. 
Resonances of the corresponding linear time-varying (LTV) system associated with the 
parametric and external excitations are identified using the method of multiple scales. 
Theoretical results are validated by available experimental results [Chapter IV]. 

1.3. DEVELOPMENT OF LINEAR TIME-INVARIANT MODELS 
1.3.1. Literature Review 

The study of geared rotor dynamics requires that the coupling between torsional 
and transverse vibrations be included in the model. Although several modeling and 
solution techniques such as lumped mass models and the use of the transfer matrix 
method have been applied to rotor dynamic problems, the finite element method (FEM) 
seems to be a highly efficient and accurate method for linear modeling. In one of the 
early examples of FEM applied to a single rotor. Nelson and McVaugh [1] used a 
Rayleigh beam finite element including the effects of translational and rotary inertia, 
gyroscopic moments and axial load. Zorzi and Nelson [2] extended this by including 
internal damping. Later, Nelson [3] developed a Timoshenko beam by adding shear 
deformation to the Rayleigh beam theory. This model was further extended by 
Ozguven and Ozkan [4] to include all possible effects such as transverse and rotary 
inertia, gyroscopic moments, axial load, internal hysteretic and viscous damping and 
shear deformations in a single model. However, none of the rotor dynamics models 
described above can handle geared rotor systems, although they are capable of 
determining the dynamic behavior of rotors which consist of shafts supported at several 
points and carrying rigid disks at several locations. 

Gear dynamics studies, on the other hand, have usually neglected the lateral 
vibrations of the shafts and bearings, and have typically represented the system with a 
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linear torsional model. Although neglecting transverse vibrations might be a good 
approximation for systems having stiff shafts, it has been observed experimentally [5] 
that the dynamic coupling between the transverse and torsional vibrations due to the 
gear mesh affects the system behavior considerably when the shafts are compliant. 
This fact directed the attention of investigators to the inclusion of transverse vibrations 
of the shafts and the bearings in mathematical models. Lund [6] developed influence 
coefficients at each gear mesh for both torsional and lateral vibrations, obtained critical 
speeds and a forced vibration response. Hamad and Seireg [7] studied the whirling of 
geared rotor systems supported on hydrodynamic bearings; torsional vibrations were 
not considered in this model and the shaft was assumed to be rigid. Iida, et al. [8] 
considered the same problem by taking one of the shafts to be rigid and neglecting the 
compliance of the gear mesh, and obtained a three degree of freedom model to 
determine the first three vibration modes and the forced vibration response due to the 
unbalance and the geometric eccentricity of one of the gears. Later, Iida, et al. 
[9,10,1 1] applied their model to a larger system which consists of three shafts coupled 
by two gear meshes. Hagiwara, Iida, and Kikuchi [12] developed a simple model that 
included the transverse flexibilities of the shafts by using discrete stiffness values, and 
studied the forced response of geared shafts due to unbalances and runout errors. They 
included the damping and compliances of the journal bearings and assumed a constant 
mesh stiffness. Although most of these gear dynamics studies have discussed several 
aspects of the problem, none of them has been able to represent a geared rotor-bearing 
system completely since almost all of them have used one or more simplifying 
assumptions such as rigid shafts, rigid bearings, rigid gear mesh etc. which may not be 
applicable to a real system. In addition, these studies proposed lower order lumped 
mass linear models. 
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Neriya, et al. [13] employed a dynamic finite element model which eliminates 
many of the simplifying assumptions. They found the forced vibration response 
of the system at the shaft frequency, excited by mass unbalances and runout errors of 
the gears by using the modal summation. But they did not consider the high 
frequency, internal, static transmission error excitation which has the major role in 
noise generation. An extensive survey of linear mathematical models used in gear 
dynamics analyses is given in a recent paper by Ozguven and Houser [14]. 

1.3.2. Mathematical Model 

Here, we assume linear bearings and no tooth separations with time-invariant 
mesh stiffness, i.e. [k] * [K(t)] and {f(q(t))} = {q(t)}. Then the corresponding LTI 

form of equation ( 1 . 1 ) is: 

[M]{q"(f)} + [C]{q'(t» +[K]{q(t)} = {F(i>} ( , 2) 

Since many investigators have modeled the generic system shown in Figure 1.1a as a 
single degree of freedom model, our analysis uses it as a reference model to transform 
the governing equation into the dimensionless form. Focusing only on the gear pair as 
shown in Figure 1.2, the equation of motion of the semi-definite system is given in 
terms of the relative translational displacement u(t) = (d gl 0g 1 -dg2eg 2 )/2 

d 2 u(t) du(t) 

™cl ~ ^ + c h — ^ + k hu(0 = F ( t ) (1.3a) 

dt dt 

where m^ is the equivalent gear pair mass defined as 
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m cl = 


1 


(1.3b, c) 





First, we establish the dimensionless time t as t = (o n t. Second, we use the base circle 
diameter of the pinion d gl and the equivalent mass m c j as characteristic length and 
mass parameters, respectively, to obtain the governing equations in the dimensionless 
form as 


[M]{q(t)} + [C]{q(t)} + [K]{q(t)} = {F m } + {F^t)} + (F aT (t)} . (1.4) 


where an overdot means derivative with respect to time t, and the dimensionless forcing 
vector consists of a mean force vector { F m } and two time-varying components: a) high 
frequency excitation due to the kinematic gear transmission error {F^ft)} and, b) other 
excitations due to mass unbalance U^, geometric eccentricities (run-out errors) Egj and 
prime mover and load torque pulsations T gia (t), typically at low frequencies, are 
combined into a single term (F aT (t)}. 

1.3.3. Gear Mesh Formulation 

The gear mesh is represented by a pair of rigid disks connected by a translational, 
viscously damped spring along the pressure line which is tangent to the base circles of 
the gears as shown in Figure 1.3. By choosing the Y axis on the pressure line and the 
X axis perpendicular to the pressure line, the transverse vibrations in the X direction 
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are uncoupled from both the torsional vibrations and transverse vibrations in the Y 
direction. The dynamic mesh force F faY at the mesh point in the Y direction can be 

written in terms of the symbols given in Figure 1 .3 from the rigid body dynamics 


fhY(t) = Ch(y g ] +- ~6 t i +Egj£2 s |COS0 t i — y g2 — ^-®t 2 

-E g2 ft S 2 cos0 t2 -eN gl Q sl cos(N gl 0 tl )) + k h (y gl + -^-0 tl 

+£ gl sin 0,i -y g 2 -Eg2sin0, 2 -ejsinCNgiOn)) (1.5) 


Here, only the fundamental component of the static transmission error e(t) is 
considered, i.e. e(t) = esin(N g iQ sl t) where N gl is the number of teeth of gear gl. 

Dynamic mesh force F faY also introduces moments T hgl and T hg2 about 
instantaneous centers of the gears: 


T hgi(*) = ^Y(O^-y- + e g i c °s0ti j 


fd, 

V 2 


T fag 2 (t) = -^iy(0( + e g 2 cos0 t2 


(1.6a) 


(1.6b) 


The total angular rotation 0 t j of the i-th gear is 0 ti (t) = fl si t + 0gj(t) , i=l,2. The 
displacement vector can be decomposed into two parts as {q} = {qs} + {qh} where 
{qs} = [*l»yi.9l. — ,x j’yj’®j» — > x n>yn>®n]’ j = lh>n, (excluding the 

degrees of freedom of the nodes where the gear and pinion are mounted) and 
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{%H x gl - Ygl >®gl * x g2 * Yg2 >®g2 ] , i=l,2. The coupling effect due to gear mesh is 
only seen in the terms governed by {%} and the mesh stiffness [K^], which 
represents this dynamic coupling due to gear mesh is obtained from equations (1.5) and 
(1.6) as 


[Kb] = k, 


0 

0 

0 

0 

0 

0 


0 

1 

fsl 

2 

0 

-1 

d g2 


0 

£gi 

2 

j2 

d gl 

4 

0 

2 

d gl d g2 


0 

0 

0 

0 

0 

0 


0 

-1 

2 

0 

1 

d g2 


0 


*g2 


d„,d 


gl_g2 


0 

}g2 

2 

j 2 

j g2 


(1.7a) 


Since mesh viscous damping is proportional to the mesh stiffness, 


N=Mk„]. 


(1.7b) 


In the overall force vector, non-zero terms correspond to displacements {qh} and are 
given in the following form: 


13 


{p*(*)}+{f.t(<)}= 


2T, 


Ugift^sinOsit + F^ 

U g Ai cosfljit 

~ E gl COSn sl t + + T gla (t) 


»*1 


U g2 Qs2 sinft S 2 f -^(t) 

u g 2^s2 cosa s2 t 


2T, 


glm 


e g 2«?2cosn s2 t-^-F 1 (t)-Tg 2a (t) 

°gl L 


(1.8a) 


Fi(t) = c h (e g2 fl s2 cosft s2 t - e gl fl sl cosQ sl t + eN gl n sl cos(N gl ft sl t)) 

+k h (e g2 sinn s2 t-e gl sinn sl t + csin(N gl a sl t)) (1.8b) 


where T g j m is the mean input torque. 

1.3.4. Finite Element Formulation and Eigen-Value Problem 

The finite element method has been employed in obtaining stiffness and mass 
matrices of the system of equations (1.4). The shafts are discretized and five degrees of 
freedom are defined at each node, only the axial motion being excluded. The stiffness 
and mass matrices of each finite rotor element are derived by using the variational 
principle [3,4,15]. The system overall matrices are obtained by combining element 
matrices, dynamic coupling matrices due to gear mesh defined by equation (1.7), 
assumed bearing damping values and stiffnesses, and the inertias of the gears and other 
lumped inertia elements. Here, [M] is diagonal and positive definite and [K] is 
symmetric and positive semi-definite. The equation describing the undamped free 
vibration of the system is obtained from equation (1.4) as 
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[M]{q(t)} + [K]{q(t)} = {0}. 


(1.9) 


The assumed form of the solution {q(t)} = {q a }cos(©t +$) is substituted into equation 
(1.9) to obtain the standard eigen- value problem: 

[K]{q a } = ® 2 [M]{q a }. (110) 

The solution of equation (1.10) yield the eigen-vectors or modes |V|/ r } and associated 
eigen-values or natural frequencies to r . Here, a torsional rigid body mode at (0=0 
exists since [K] is semi-definite. 

1.3.5, Forced Vibration Response 

The excitation given by equation (1.8) (for no torque pulsations, i.e. 
T gla (t) = Tg 2a (t) = 0) is the sum of three sinusoidal terms at frequencies ft s j, fi S 2 and 

gear mesh frequency NgjQ s i 

{F(t)} = X {Fgj } sin^gjt + <t> si ) + {1^ }sin(N gl O sl t + 4>h) (1 ' 1 1 ) 


The steady state displacement response of the system due to this excitation is assumed 
to be 


{q(t)} = X [p gi ]{F gi }sin(n si t + (J» si ) + [P h ]fa} sin ( N gAl' + <t>h) (112) 

i=l L 
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where [Pgj] i=l, 2 and [Pt,] are the dynamic compliance matrices in the frequency 
domain corresponding to the excitation frequencies, ft sl ,ft s2 and N gl ft sl , 
respectively. 


f8 1= y (VrHVrr 


i = 1. 2; 


(1.13a) 


IM- 1 


r=l 


{^rHVr) 


®r (Ngifigi) ^ + 2j^ r Ngif2 si a) r 


(1.13b) 


I Vr) represents the r-th mass matrix -normalized modal vector, n is the total number of 
degrees of freedom of the system, j = V-l, and £ r is the damping ratio for r-th mode. 
See Appendix A for the computer code GRD which uses the theory given here. 

1.4. PARAMETRIC STUDIES 
1.4.1. Modes of Interest 

The generic system shown in Figure 1.1 is modeled by FEM to examine the 
natural modes of a general linear geared rotor-bearing system and to study the effects of 
several system parameters on the dynamics of the system. Two numerical data sets as 
listed in Table 1.1 are used. Predicted natural frequencies and modes for Set A are 
presented in Table 1.2. The response of the system to e(t) is also computed; Figures 

1.4, 1.5 and 1.6 display the response in the Y and torsional directions at the pinion 
location and the dynamic load to static load ratio at the mesh point d gl ^j Y / 2 T mg i, 

respectively. The system has no peak responses at the modes corresponding to motion 
in the X direction since the excitation is applied in the Y (pressure line) direction, and 
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Table 1.1 Numerical data sets of the system used for calculations. 


Parameters 

Set At 

Set B 

I g l,I g 2 (kg-m 2 ) 

0.0018, 0.0018 

0.0097, 0.0097 

m gl’ m g2 ( k g) 

1.84, 1.84 

3.45, 3.45 

d g i, d g 2 (m) 

0.089, 0.089 

0.135, 0.135 

N gl 

28 

30 

k h 

1.0x10 s 

1.0x10 s 

k b (N/m) 

variable 

rigid 

kslL’ L slR ( m ) 

0.127, 0.127 

variable 

L s2L’ L s 2R ( m ) 

0.127, 0.127 

variable 

K 1 ,K 2 (N-m/rad) 

— 

variable 

I P , I L (kg-m 2 ) 

— 

variable 

^sl,0’ ^s2,o ( m ) 

0.037, 0.037 

0.04, 0.04 

^sl,i’ ^s2,i ( m ) 

0., 0.01 

— 

e (m) 

9.3xl0' 6 

— 


t NASA Lewis Research Center gear test rig. 
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Table 1.2. First 10 natural frequencies for set A of Table 1.1 (Jq/k^lO). 


Modal 

Index 

Modes of 
Interest 

Natural Frequency o\ Natural Mode 

Hz. Description 

0 


0 

torsional rigid body 

1 

Vi 

581 

first transverse-torsional coupled 

2 


687 

X direction, transverse (driving shaft) 

3 

Vh 

689 

Y direction, transverse 

4 


691 

X direction, transverse (driven shaft) 

5 

Vin 

2524 

second transverse-torsional coupled 

6 


3387 

Y direction, transverse 

7 


3387 

X direction, transverse (driving shaft) 

8 


3421 

X direction, transverse (driven shaft) 

9 


3421 

Y direction, transverse 
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response in Y -direction (m) 



excitation frequency (Hr) 


Figure 1.4. Forced response of the system with dataset A (Table 1.1) to the 
displacement excitation in the direction of pressure line (at pinion 
location) for four different bearing stiffnesses k^. 
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excitation frequency (Hz) 


Figure 1.5. Forced torsional response of the system with dataset A (Table 1.1) to 
the displacement excitation at pinion location for four different bearing 
stiffnesses k b . 
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dynamic to static load ratio 



Figure 1.6 Dynamic to static load ratio for dataset A (Table 1.1) due to the static 
transmission error excitation for four different bearing stiffnesses kj,. 
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the vibration in the X direction is dynamically uncoupled from the vibrations in the Y 
and torsional motions, as described in the previous section. Therefore, the natural 
modes corresponding to motions in the X direction can be eliminated when only e(t) 

excites the system. Accordingly, the following three typical modes are of special 
importance: 

I. First Transverse-Torsional Coupled Mode yp This is the second mode listed 
in Table 1.2 at C0j (581 Hz for the system considered). This mode corresponds to the 
first peak in Figures 1.4 and 1.5 and its schematic shape is shown in Figure 1.7a. 
Here, shafts move in opposite directions, gears vibrate in opposite directions also; but 
transverse and torsional vibrations combine to yield small relative motion at the gear 
mesh point. Therefore, dynamic loads at the mesh point are not large, resulting in no 
peak in Figure 1 .6 at o>] while Figures 1 .4 and 1.5 have peaks governed by this mode. 

II. Purely Transverse Mode At this mode with natural frequency o>jj (689 
Hz in Table 1.2), there is no torsional vibration and both shafts vibrate in phase in the 
pressure line direction as shown in Figure 1.7b schematically. The relative 
displacement at mesh point is zero since the gear ratio v g =l. Therefore, e(t) cannot 

excite this mode as no peaks are observed in Figures 1.4, 1.5 and 1.6 at frequency (Dy. 

III. Second Transverse-Torsional Coupled Mode Ym : The second and the 
highest peak seen in Figures 1.4, 1.5 and 1.6 corresponds to this mode at (Dm . As 
seen from the mode shape illustrated in Figure 1.7c, both shafts and gears vibrate in 
opposite directions, and transverse and torsional vibrations are additive at the mesh 
point. Thus a large relative displacement at the mesh point is obtained, which results in 
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Figure 1.7. Typical natural modes of interest; see Table 1.2 for further details.. 
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a large peak in Figures 1.4, 1.5 and 1.6. This is the mode at which the coupling 
between transverse and torsional vibrations is very strong. 

These three modes are observed in all geared rotor systems and they play an 
important role in governing dynamic response of the system excited by e(t). 

Therefore, any accurate mathematical model of the geared rotor-bearing systems must 
be able to predict these modes. 

1.4.2. Effect of Bearing Compliances 

In most cases, radial stiffness of a typical rolling element bearings k b is roughly 
of the same order of magnitude as the gear mesh stiffness k,,; in general it lies in the 
range 0.1k h <k b <100k h . Therefore, bearing flexibility should be included in the 

analysis. A parametric study for data set A of Table 1.1 has been conducted to 
demonstrate the effect of k b on the natural frequencies and the frequency response of 
the system excited by e(t). Figures 1.4, 1.5 and 1.6 show the response in the Y and 
torsional directions and the dynamic mesh load respectively for kh values ranging from 
0.1k h to 10k h , shaft stiffnesses are being kept the same. An increase in bearing 

stiffness results in an increase in both natural frequencies and the peak amplitudes. 
Above the range of kj,, the bearing becomes very rigid when compared to the shaft 

compliances and its effect on the system can be ignored. 

1.4.3. Effect of Shaft Compliance 

Data set B of Table 1 . 1 has been used to study the effect of shaft compliance on 
the natural modes. The shaft length is varied and the corresponding co r values are 
predicted by FEM; Figure 1 .8 displays this result. The first natural frequency Oj does 
not change considerably with varying shaft length, whereas ©u and (Ojn are strongly 
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dimensionless natural frequency 




dependent on the shaft length, especially at smaller lengths. Another observation from 
Figure 1 .8 is that cojj and ©m become very large for smaller shaft lengths and clearly 

move beyond the range of the operational speed of most geared rotor systems. In this 
case, the first coupled transverse-torsional mode Vj can be assumed to be uncoupled 
from these two modes. 

1.4.4. Effect of Load and Prime Mover Rotary Inertias 

As shown in Figure 1 1 , the shafts are connected to a prime mover and a load at 
either end through flexible torsional couplings. The following parameters need to be 
considered: motor and load rotary inertias, torsional compliances of the flexible 
couplings, and stiffnesses of the driving and driven shafts. 

First, the motor and load are assumed to be connected to the shafts without 
considering any torsional couplings in between. Figure 1.9 displays the variation in 
311(1 with a variation in the prime mover inertia. Here, data set B of Table 
1.1 with L si =0.04 m is used, ©n and ©m are not affected and therefore, the inertias of 

motor and load can be disregarded if the major concern is to predict these two modes. 
However, ©j is strongly dependent on prime mover and load inertias. 

Second, the load and prime mover inertias are fixed (I p =5I gl ) and the torsional 
springs Kj and K 2 , which represent flexible couplings and shafts are varied. Figure 
1.10 shows the variation in ©, with changing K, and K 2 . Here ©u and ©m are again 
not affected, as expected. On the other hand, ©j is almost constant (which is nearly 
equal to the value yielded by zero prime mover and load inertias) up to a point and, it 
starts increasing with increasing K. This indicates that the motor and load are isolated 
from the geared rotor system when the in-between torsional elements are compliant 
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dimensionless natural frequency 



2 : 



dimensionless natural frequency 



torsional stiffness (N.m/rad) 


Figure 1.10. Effect of the torsional stiffness K of the transmission elements in 
between the gear box and the prime mover & load inertias on the natural 
frequencies. 
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enough, which is the case in most practical systems. Under these circumstances, 
motor and load inertias can be neglected in the analysis. 

1.5. REDUCED ORDER LINEAR TIME-INVARIANT MODELS 

In this section, three different reduced order analytical models of the geared rotor- 
bearing systems shown in Figure 1.1 will be developed and the conditions and system 
parameters at which these simple models can predict the dynamics of the system 
accurately, will be discussed. The finite element model of Section 1.3 will be 
employed as a reference model to check the validity. 

1.5.1. Single Degree of Freedom Torsional Model of Gear Pair 

As the simplest model, a single degree of freedom (SDOF) model of geared rotor 
systems shown in Figure 1.2 is considered. The shaft and bearing flexibilities and the 
motor and load inertias are not considered in this model. This model can only predict a 
single mode at a frequency (0 n as defined in equation (1.3c) which corresponds to the 
first transverse-torsional mode at CDj . Here, 0 >j— ><o n when shaft lengths L s j— >0 and 
bearings are very stiff. And, C 0 j] and cojjj are sufficiently beyond the operational speed 
range and the variation in (Oi is assumed to be negligible with shaft length as shown in 
Figure 1.8. Therefore, in some cases, the SDOF model of Figure 1.2 can be used to 
represent the system provided these conditions are met. For instance, for data set B of 
Table 1.1 with 1^=5 cm, a SDOF model can be utilized up to an operational rotational 

speed of 6000 rpm which corresponds to the excitation frequency at 3000 Hz for 
N gl =30 teeth. As it is seen in Figure 1.8, only the first mode is observed and its 
variation is not significant within 0^0<3000 Hz and O^L^^S cm. 
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1.5.2. Three Degree of Freedom Model 

The SDOF model of Figure 1 .2 is not adequate when the shafts and bearings are 
compliant. To overcome this deficiency, a three degree of freedom (3-DOF) 
transverse-torsional model as shown in Figure 1.11a is developed. Equation 1.4 gives 
the equations of motion with dimensionless mass [M], damping [C] and stiffness [K] 
matrices and displacement vector (q) given as follows: 


[M] = 


1 

0 

0 


0 

rrigi / m cl 
0 


0 

0 

m g2 / m cl 


(1.14a) 


[C] = 


Cb 


V k h m cl 


1 

1 

-1 


1 -1 
(1 + C(j S | / Cjj ) —1 

_1 ( 1 + Cbs2/ C b). 


(1.14b) 


[K] = 


1 

1 

-1 


1 -1 

( 1 + k bsl/ k h) “1 

-1 (l + k bs2 /k h ) 


(1.14c) 


{q}=[ u - y g i> y g 2] 


(1.14d) 


u _9 g i d g2 9 g2 

2 2d gl ’ 



i = l,2. 


(1.14e,f) 


Here, k bsl and k bs2 are equivalent lateral stiffnesses representing shaft and bearing 
flexibilities, and c bsl and c bs2 are equivalent viscous damping values. This model can 
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Pinion 



Figure 1.11. Reduced order analytical models of Figure 1.1; a) three degree of 
freedom model, b) six degree of freedom model. 
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accurately predict all three modes of interest as evident from Figure 2.8 and Table 1.3 
where its predictions are compared with the results of FEM. 

When the system is connected to the motor and load inertias, then the in between 
torsional stiffnesses Ki and K 2 should be compliant enough to be able to neglect the 
effects of the motor and load inertias, as it is discussed earlier in Section 1.4.4. In 
summary, the 3-DOF model shown in Figure 1.11a can be used to describe the 
dynamics of the geared rotor system when: a) the shafts and bearings are compliant and 
provided the shafts are short such that higher order bending modes of the shafts are out 
of the frequency range considered, and b) the torsional stiffnesses of the connections in 
between the motor and load inertias and the gear box are sufficiently compliant. 


1.5.3. Six Degree of Freedom Model 

A six degree of freedom (6-DOF) model as shown in Figure 1.11b can be 
employed to represent the geared rotor-bearing system when the effects of the motor 
and load are not negligible as mentioned in Section 1.4.4. Equations of motion are still 
given by equation (1.4) and the dimensionless system matrices are defined as 


[M] = 


m cl 


Ip/d|x 

0 

0 

0 

0 

0 


0 

[ g ,/dii 

0 

0 

0 

0 


0 

0 

Ig2/ d «l 

0 

0 

0 


0 

0 

0 

I L /d|l 

0 

0 


0 

0 

0 

0 

m gl 

0 


0 
0 
0 
0 
0 

m g 2j 


(1.15a) 
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Table 1.3. Comparison of typical modes obtained by FEM and 3-DOF models for set 
B of Table 1.1; ^^,=0.3. 


Mode \\f 


Vi 

Vu 

Vrn 

Displacement 

FEM/3-DOF 

FEM/3-DOF 

FEM/3-DOF 

Ggl 

-0.012/-0.011 

-0.001/0.0 

-0.210/-0.246 

e g2 

0.012/0.011 

-0.001/0.0 

0.210/0.246 

y g i 

1.0/1 .0 

1.0/1 .0 

-1.0/- 1.0 

yg2 

-1.0/- 1.0 

1. 0/1.0 

1.0/1 .0 
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[C] = 


V k h m 


cl 


IK1 ^ 


Cj_ 

d gl 


Kj_ 

d gi 


°gi 

{ ^L + £h) 

‘oil 4 


0 

c h v g 

4 

(% + ^ 

d gi 4 


) - 


-symmetric - 


-symmetric - 


C 2 


C2_ 

d gl 


Ki 

'4. 

0 

0 

<^> 
»!> 4 

k h v g 

0 

4 




k 2 


k 2 


J gi 


0 


0 — 


c h 
2 

c h v g 

2 

0 

(Ch +C b.l) 


0 

c h 
2 

c h v g 
2 

0 

"Ch 
(Ch+Cbs2>J 


(1.15b) 


0 

^h_ 

2 

k h v g 

2 

0 

(kh + kb,i ) 


0 

_^h_ 

2 

Vg 

2 

0 

~ k h 

(k h +k bs2 )J 


(1.15c) 


{q} = [®p- e gi* 0 g2. e L> ygi* y g 2] 


(1.15d) 


For data set B of Table 1.1 predictions yielded by the 6-DOF model are compared 
with those by FEM as shown in Figure 2.9 and Table 1.4. Based on these results, it 
can be concluded that the 6-DOF model is accurate enough to predict the natural modes 
of the system. Accordingly, 6-DOF model must be employed when the effects of the 
motor and load are not negligible. 
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Table 1.4. Comparison of typical modes obtained by FEM and 6-DOF models for set 
B of Table 1.1; Lsj/d gl =0.3, Ip/I g i=5. 


Displacement 


0p 

0 gl 

0 g2 

ygi 

yg2 


Mode tjr 


Vi 

FEM/6-DOF 

Vn 

FEM/6-DOF 

Vm 

FEM/6-DOF 

0.079/0.079 

0.002/0.0 

-0.018/-0.013 

-1.0/- 1.0 

-0.005/0.0 

1.0/1 .0 

1.0/1 .0 

-0.005/0.0 

-1.0/- 1.0 

-0.079/-0.079 

-0.002/0.0 

0.018/0.013 

0.006/0.006 

1.0/1 .0 

0.391/0.451 

-0.006/-0.006 

1.0/1 .0 

-0.39 1/-0 .451 
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1.6. CONCLUSION 

In this chapter, a finite element model to investigate the dynamic behavior of 
linear time-invariant geared rotor systems has been developed. The transverse 
vibration of the system associated with shaft and bearing flexibilities and the dynamic 
coupling between the transverse and torsional vibrations due to gear mesh have been 
considered. Natural modes of the system have been identified and forced vibration 
response due to both low frequency external and high frequency internal excitations 
have been determined. Reduced order analytical models of the geared rotor bearing 
system have also been developed. Three different linear time-invariant models (SDOF, 
3-DOF and 6-DOF) have been suggested to represent the geared system. By 
comparing results with FEM predictions, it has been shown that such reduced order 
linear models are reasonably accurate. Therefore these models will be extended in 
Chapters II, IQ and IV in analyzing the effects of system non-linearities and time- 
varying mesh stiffness. 
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CHAPTER II 

NON-LINEAR DYNAMIC ANALYSIS OF 
A SPUR GEAR PAIR 


2.1. INTRODUCTION 

2.1.1. Excitation Types and Backlash 

The focus of this chapter is on the backlash non-linearity as excited primarily by 
the transmission error between the spur gear pair. A gear pair is bound to have some 
backlash which may be designed to provide adequate lubrication and eliminate 
interference due to manufacturing errors. Backlash-induced torsional vibrations may 
cause tooth separation and impacts in unloaded or lightly loaded geared drives. Such 
impacts result in intense vibration and noise problems and large dynamic loads, 
which may affect reliability and life of the gear drive [16,17]. Excitation mechanisms 
can be grouped as follows: 

A. External Excitations: This group includes excitations due to rotating mass 
unbalances, geometric eccentricities, and prime mover and/or load torque fluctuations 
[18]. Although mass unbalances and geometric eccentricities can be reduced through 
improved design and manufacturing, torque fluctuations are not easy to eliminate 
since they are determined by the characteristics of the prime mover (piston engines, 
dc motors etc.) and load [19]. Such excitations are typically at low frequencies 

which are the first few multiples of the input shaft speed Q s . Practical examples 
include rattle problems in lightly loaded automotive transmissions and machine tools 
[19,20], 
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B. Internal excitations; This group includes high frequency excitations 

caused by the manufacturing related profile and spacing errors, and the elastic 
deformation of teeth, shafts and bearings. Under the static conditions, all such 
mechanisms can be combined to yield an overall kinematic error function known as 
"the static transmission error" e(F) [17,18]. This error is defined as the difference 
between the actual angular position of the driven gear and where it would be if the 
gears were perfectly conjugate [17,18,21-23]. In gear dynamic models, e(F) is 
modeled as a periodic displacement excitation at the mesh point along the line of 
action [15,24-26] and its period is given by the fundamental meshing frequency 
Q fa =N Q s where N is the number of teeth on the pinion. Practical examples 

include steady state noise and vibration problems in automotive, aerospace, 
industrial, marine and appliance geared systems. 

2.1.2. Literature Review 

Experimental studies on the dynamic behavior of a spur gear pair with backlash 
started almost 30 years ago and still continue [27-29]. As one of the better examples 
of such experiments, Munro [27] developed a lightly damped (damping ratio £=0.02) 
four-square test rig to measure the dynamic transmission error of a spur gear pair. 
He used high precision gears with rigid shafts and bearings, and showed 
experimentally that the tooth separation takes place when the mean load is less than 
the design load. Dynamic transmission error versus speed curves were plotted to 
illustrate the steady state response and the jump phenomenon. Kubo [28] measured 
the dynamic tooth stresses using a similar set-up in order to calculate the dynamic 
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factors. He also observed a jump in the frequency response of the gear pair with 
backlash even though the test set-up was heavily damped (£=0.1). 

Such experimental studies, though limited in scope, have clearly shown that the 
gear pair dynamics can not be predicted with a linear model - see Ozguven and 
Houser [14] for a detailed review of the linear gear dynamic models as available in the 
literature. Although most of the non-linear mathematical models used to describe the 
dynamic behavior of a gear pair are somewhat similar to each other, they differ in 
terms of the excitation mechanisms considered and the solution technique used. For 
instance, a large number of studies have focused on the rattle problem in lightly 
loaded geared drives which are excited by the low frequency external torque 
excitations [30-35]. A few investigators have included the static transmission error 
excitation in the non-linear models [24,35-37]. 

The gear backlash non-linearity is essentially a discontinuous and non- 
differentiable function and it represents a strong non-linear interaction in the 
governing differential equation. This issue has been discussed by Comparin and 
Singh [33] and they have concluded that most of the solution techniques available in 
the literature can not be directly applied to examine this problem. Most of the gear 
dynamic researchers recognized this problem implicitly, and therefore employed 
either digital or analog simulation techniques [19,24,29,35-38]. For instance, 
Umezawa et.al [29], Yang and Lin [32] and Ozguven and Houser [24] have solved a 
one degree of freedom torsional model of the gear pair using numerical techniques. 
Lin et.al. [38] included motor and load inertias in a three degree of freedom torsional 
model. Kucukay [35] has developed an eight degrees of freedom model to include 
the rocking and axial motions of the rigid shafts. In most of these studies, with the 
exception of Umezawa's analysis [29] which did not include any backlash, a 
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discontinuity has been seen in the frequency response characteristics. But many 
investigators have typically joined two discrete points to show a broad jump in the 
frequency response curve [24,35,38]. Some of these problems have been due to the 
numerical simulation techniques which may not work or may result in misleading 
answers if not employed properly. Such difficulties have been found by Comparin 
and Singh [33], Singh et al. [34] and Gear [39-40] but are yet to be resolved or 
addressed by the gear dynamics researchers. Accordingly, one of the major 
objectives of this chapter is to examine whether numerical simulation techniques can, 
in fact, be used to predict the dynamic response completely, and what precautions one 
must take to develop such a mathematical model. Since Comparin and Singh [33] 
and Singh et al. [34] have examined the external excitation problem, in this chapter 
we focus mainly on the internal excitation and see whether the numerical simulation 
technique can be made to work for the prediction of the non-linear frequency 
response characteristics. 

A few researchers have attempted to obtain the analytical solutions for a gear 
pair problem, based on the piece- wise linear techniques which divided the non-linear 
regime into several linear regimes [41-43]. For instance, Wang [41,42] has used two 
and three degree of freedom torsional models with backlash, and assumed that the 
gear teeth are rigid and the driven gear has an infinite inertia. The governing 
equations have been solved using the piece-wise linear technique. It should be noted 
that the piece-wise linear technique gives only solutions for the equivalent linear 
systems and one typically may have difficulties in combining such solutions [43]. 
Comparin and Singh [33] overcame these problems by employing the harmonic 
balance method (HBM) and constructed analytical solutions for the non-linear 
frequency response characteristics of a gear pair with backlash as excited by the 
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external torque. In this chapter, we will use the same technique to examine the 
internal excitation problem, and compare results with digital simulation and 
experimental studies. Further literature review is included in subsequent sections. 


2.2. PROBLEM FORMULATION 
2.2.1. Physical Model 

A two degree of freedom semi-definite model of the spur gear pair with rotary 
inertias I gl and I g2 and base circle diameters d gl and d g2 as shown in Figure 2.1 is 

considered here. The shafts and bearings are assumed to be rigid. The gear mesh is 
described by backlash of 2b and by a time invariant mesh stiffness kjj * kj,(t) 
when in contact and viscous damping c fa . The equations of torsional motion of the 

gear pair as shown in Figure 2.1 are 
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(d . 
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Figure 2. 1 . Gear pair model. 
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where T tl (r)=f glm + T gU (D end tJD = T g2m + f f JT) are .orques 

on pinion and gear and f is a non-analytical function essentially describing the mesh 
elastic force as shown in Figure 2.1. Here, output torque fluctuation T g ^(t) will 

be neglected to simplify the dynamic problem, i.e. T g2 (^) Equations 

(2.1a) and (2.1b) can be reduced to one equation in terms of q(D which is defined 
as the difference between the dynamic transmission error x(t") and the static 
transmission error e (F). 

m c ,^ + c 11 § + k b f(q(i)) = F n , +F aT (i)-m c ^f ; (2.2a) 

dt dt at 

5(F) = x(F) - e(F) = %e g| (D - ^-e g2 (D - 1(F) ; (2.2b) 


= ( d 2 


gl 

K. 
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g2 ; 
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g_L . 

m c2~ d ’ 
gl 


(2.2c, d) 


2T 


F m = 


glm 


2T 


g2m 

'«* 


2m c T - (F) 

V«> = STX 


c2 g 1 


(2.2e,f) 


f(q(D) = 


f(q(0 ) 



q(F) -b; 

^ 0 ; 

q(t ) + b; 


q(F) > b 
- b < q(D < b • 
q(D < - b 


(2.2g) 
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where m cl is the equivalent gear pair mass, F m is the average force transmitted 
through the gear pair, F aT (t ) is the fluctuating force related to the external input 

torque excitation and f(q(t )) is the nonlinear displacement function. Equation 
(2.2a) is nondimensionalized by letting q(F) = q(F)/b, © n = /m d t 
t = (o n t and C= c h / ( 2m c i 0) n) • Now * consider harmonic excitation for both 
e(D and F aT (F) as e(F) = e sin(Q h F + <|> h ), F aT (F) = F #T sin ( H y F +<|> T ) 
where and 12^, are the fundamental excitation frequencies of internal 

displacement and external torque fluctuations, respectively. Further, define 
dimensionless excitation frequencies f2 h = f2 h /to n and Q T = f2 T /co n , 

dimensionless external mean load F m = F m /bk h , amplitudes of the dimensionless 

internal (F^t)) and external alternating forces (F_(t)) F.=F,,/bk L and 

c aT aT h 

F ah = * ^ b and nonlinear displacement function f(q(t)) to yield the following 
governing equation of motion. 


q(t) + 2£q(t) + f(q(t)) = F(t) 


(2.3a) 


F (t)= F m + F aT si n (Q T t +V + F * ft h sin (*V + 4> h ) (2.3b) 


f(q(t)) = 


f(q(t)) 

b 


q(t) - 1; q(t) > 1 
- 0 ; - 1 < q(t) < 1 

q(t) + 1; q(t) < — 1 


(2.3c) 
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2.2.2. Scope and Objectives 

When only external forces excite the system, i.e. Fj(t)=0, equation (2.3a) 
reduces to 

q(t) + 2 ^q(t) + f(q(t)) = F m + F e (t) = F m + F #T sin (Q T t + <t» T ) (2.4) 

This equation has been solved both analytically and numerically by Comparin and 
Singh [33]. Conversely, no analytical solution is available when the system is 
excited by internal static transmission error at the mesh frequency ft h which is 
considerably higher than O t . The governing equation is given by substituting 
F e (t)=0 in equation (2.3a); note that the external mean load F m is not equal to zero. 

q(t) + 2Cq(t) + f(q(l)> = F m + F .(t) = F m + sin (fl h t + * b ) (2.5) 

Both equations (2.4) and (2.5) include the clearance non-linearity. While 
equation (2.4) represents the conventional representation of the vibro-impact problem 
[33,34], equation (2.5) is more applicable to the clearance problems in built-up 
assemblies where the excitation is generated by the kinematic errors. This chapter 
focuses on the steady frequency response characteristics of equation (2.5) which 
represents a gear pair with backlash as excited harmonically by the static transmission 
error excitation e(t ) or F^. Specific objectives of this chapter are as follows: 

1. Solve equation (2.5) numerically to resolve various modeling issues such as 
the existence of multiple solutions, subharmonic resonances and chaos. 
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2. Construct analytical solutions to equation (2.5) using the harmonic balance 
method (HBM) which has been applied successfully to solve equation (2.4) by 
Comparin and Singh [33J. 

3. Compare digital simulation and harmonic balance techniques and establish 
the premises under which the jump phenomenon can be predicted. 

4. Perform parametric studies in order to understand the effects of F ,F. and 

nv ab 

C the frequency response. Vary the force ratio F = F m /F a j J which is a measure of 

the load on the gear pair and compare the dynamic behavior for lightly and heavily 
loaded gears. 


5. Validate analytical and numerical solution techniques by comparing these 


with previous experimental studies [27,28]. 

6. Compare the frequency response characteristics of equations (2.4) and (2.5), 
and also examine the possibility of finding overall response when both external and 
internal excitations are applied simultaneously. 

7. Consider the periodic static transmission error excitation case, i.e. 


F(t)= F m + F.(t)= F m + £ (ja h ) 2 F ah .sin(jn h t + <t» hJ ); only 


the first 


three (k=3) harmonics are included. 


2.3. DIGITAL SIMULATION 

Clearance or vibro-impact problems in single degree of freedom systems have 
been examined by a number of investigators whose formulations are similar to 
equation (4) - see Comparin and Singh [33] for a detailed review. Moreover, Shaw 
and Holmes [43] and Moon and Shaw [44] have considered an elastic beam with one 
sided amplitude constraint subject to a periodic displacement excitation, and have 
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shown experimentally and numerically that the chaotic and subharmonic resonance 
regimes exist. Whiston [45-47] has investigated the non-linear response of a 
mechanical oscillator preloaded against a stop. He has solved die system equation for 
harmonic excitation by digital simulation and studied the existence and stability of the 
subharmonic and chaotic responses and the effect of preload on chaos. Similarly, 
Ueda [48] has solved the Duffing's equation, q + 2£q + q^ = F sin t , numerically 

and defined the regions of different solutions on a £ versus F map. According to 
him, the existence of harmonic, subharmonic and chaotic responses depends on 
values of C, and F, and multiple steady state solutions typically exist. Thompson and 
Stewart [49] have reviewed the available literature, with focus on the Duffing's 
equation. It should be noted that equation (2.5) is different from the non-linear 
differential equations considered by the above mentioned studies. Therefore, 
equation (2.5) must be studied in depth as the results of the other non-linear equations 
may not be directly applicable to our case. 

First, we solve the governing non-linear differential equation (2.5) numerically 
using a 5th-6th order, variable step Runge-Kutta numerical integration routine 
(DVERK of IMSL [50]) which is suitable for a strongly non-linear equation 
[33,39,40]. Second, we investigate the existence of chaos and subharmonic 
resonances. Since the steady state response of the system due to the sinusoidal 
excitation is of major interest, it is necessary to run the numerical program for a 
sufficient length of time. The number of cycles of the forcing function required to 
reach the steady state depends on £. 

A lightly loaded system (F=F m /F ah =0.5) with low damping (£=0.02) is 
considered as the first example case. For F m =0.1 and F a j 1 =0.2, the response q(t) is 
computed over the frequency range 0 < < 1.5. Figure 2.2 shows phase plane 
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plots q(t) vs. q(t) for different Q h values. For ft h =0.3, all transients converge to 
one periodic solution at the fundamental frequency of the forcing function 
irrespective of the initial conditions q(0) and q(0) . Therefore, it is called a "period- 
one, tp ” attractor where t p =2n/£2 h . But, in the case of Q h =0.5, three coexisting 
period-one attractors have been found as shown in Figure 2.2. Here, q(0) and q(0) 
define three steady state limit cycle solutions. For all initial conditions given by 
- 2 < q(0 ) < 2 and — 2 < q(0 ) < 2 , a map of the domains of attraction for each 
steady state solution is obtained in the map in Figure 2.3. If a smaller increment is 
used for the initial conditions, a finer resolution will be obtained. Hence each phase 
plane plot shown in Figure 2.2 is strictly governed by a subset of the initial 
conditions. Similarly, two period-one attractors are found at Q h =0.6. At ft fa =0.7, 

besides two period-one t p attractors, two more solutions of period 3t p exist, i.e. 
period-three attractors. But only a period-two 2tp attractor is seen when Q h =0.8. 
Within the range 1. 0 £ £2 h £ 1. 5 , non-periodic, steady state or chaotic response is 

observed. Figure 2.4 shows the chaotic time history and the Poincare map (strange 
attractor) at £2^= 1.0. These results are qualitatively, but not quantitatively, similar to 

the studies reported on the Duffing's equation [48,49] and clearance non-linearity 
(equation (2.4) type) [33]. 

It is concluded that the subharmonic response of period ntp provided n ^ 1 and 
the chaotic response ( t p — > «*>) are seen in the gear pair only for a certain set of 

parameters F m , F , £ and £2 h . It must also be noted that the multi-solution regions 
are strongly dependent on the choice of q(0) and q(0) . Only one steady state 
solution can be found via digital simulation when only one set of initial conditions is 
chosen at a given Q h , and the rest of the steady state solutions are not predicted. This 
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initial velocity 



Figure 2.3. Domains of attraction for three attractors; x no impact (Figure 2.2bl), • 
single sided impact (Figure 2.2b3), ° double sided impact (Figure 2.2b2); 

F =0.5, F m =0.1, C= 0.02 and fl h =0.5. 
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results in an incomplete frequency response description. These issues will be 
discussed further in Section 2.5. 

A heavily loaded system with F =2, F m =0.1 and the same amount of damping 

£=0.02 is considered as the second example case. Figure 2.5 shows the phase plane 
plots for the same values of which are used in the first example case. Unlike the 

first example case, no chaotic responses are found here. All of the solutions are 
period-one tp type for all Ojjd.O and Figure 2.6 shows the typical domains of 
attraction at O h =0.6. However, within the range 1. 0 £ Q < 1. 5 , one period-two 

attractor co-exists with the period-one orbit. Hence the force ratio F determines the 
existence of the chaotic and subharmonic responses. To illustrate this point, consider 
chaos shown in Figure 2.4. As F and £ are increased, significant changes in the 
response are observed in Figures 2.7 and 2.8. A transition from the chaos to a 
period-two, and then to a period-one steady state solution is seen when F is 
increased from 0.75 to 1.0 and then to 1.5. Similarly, an increase in £ to 0.05 
reduces chaos to a period-eight attractor, which then bifurcates to a period-two orbit 
at £=0.1 as shown in Figure 2.8. Since most real geared systems are heavily loaded 

A 

with a high F , chaotic and subharmonic responses should not be seen under the 
normal driving conditions. This issue will be discussed again in section 2.6. 

2.4. ANALYTICAL SOLUTION 

An approximate solution for equation (2.5) is constructed using the harmonic 
balance method (HBM). Assume that q(t) = q m + q a sin(Q h t + <|> r ) where q,,, and q a 

are mean and alternating components of the steady state response, and <|> r is the 

phase angle. Here, higher harmonics of the response are not included in the analysis. 
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sided impact (Figure 2.5c 1); 
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- 2.0 

2.0 

q 



A A 

Figure 2.7. Time histories of q(t) for F ah =0.2, £=0.02 and a) F=0.75, b)F = l, 
c)F = 1.5. 
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Figure 2.8. Time histories of q(t) for F m =0.1, F =0.5 and a) £=0.05, b) £= 0. 1 . 
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The quasilinear approximation to the nonlinear function f(q) with the excitation 
F(t)= F m + F (t)= F m + F^fl 2 sin (£l h t + <|> h ) is in the form [33,51] 

f(<l)= N n.<lm +N a<la sin ( ti h ,+l| ’r> + N » < l. cos ( Q l, t+ 'l l r > < 26 > 


where the describing functions N m , N a and N* are defined as: 


N 


m 


1 

2ttq m 


2 n 


|f(q m + q a sin <p> d <p ; 

0 


(2.7a) 


N 


a 


i 

*q a 


2n 


}f(q m +q a sin (p)sin <pdtp ; 
o 


(2.7b) 


N* = rco” ^ f( q*« + q» sin 9) cos< pdq> ; <P=^ h t + < t>r- (2.7c, d) 

Ha o 


Equation (2.3b) is substituted into equations (2.7a), (2.7b) and (2.7c) to obtain 


N m = 1_ T^[S(Y + )-g(Y-)] ; N a = 1- ^[h(Y + )+h(y_)] ; (2.8a,b) 


1±q «" 

N* a = 0 ; Y ± = -q— 


(2.8c,d) 


where 
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g(Y) = 


^•(y sin *Y + Vi-Y 2 ) I lYl^l. 

UyI; IyI > 1 


(2.8e) 


h(Y) = 


!'(sin- 1 Y + YVl-Y 2 ) ; lYl^l 

- 1 ; y <- 1 

1 ; Y> 1 


(2.80 


Comparin and Singh [33] used truncated series expansions for functions g(y) and 
h(y) given by equations (2.8e) and (2.80- They stated that the error involved in 
using truncated series is within 6 percent when only the first two terms are considered 
with the coefficient of the second term adjusted to yield the actual value for y = 1 . 
Using the same approach, one gets 


g(y )s n( 1+ (~2~)y 2 ) ; h (y> s ^(y~ (■^~4^')y 3 ) ; m*i 


and obtains the following frequency response by substituting equations (2.6) into 
equation (2.5) and equating the coefficients of like harmonics: 


q a = 


F M 2 
ah h 


V( n *-° b) J+ 




N, 


(2COJ' 


i r = 4> -tan 


- 1 


2CQ. 


N a-^ h y 


(2.10a,b) 


(2.10c) 
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Depending on the damping ratio £ and the parameters F m , and £2 which 
define the excitation, there are three cases at which different solutions are obtained: 
(a) no impact (no tooth separation), (b) single sided impacts (tooth separation, but no 
back collision), and (c) two sided impacts (back collision). 

Case I: No Impact : The tooth separation (impact) is not observed in a geared 
system if the displacement q(t) lies in the region q(t)>l all the time. This condition, 
shown in Figure 2.9 as case I, can be described mathematically as 

|qm + 9a| > l and I t lm -< la| >1 (211) 

Then, for no impact region in which the conditions defined by equation (2.1 1) are 
satisfied, the describing functions are given by 

N = i — • N = 1 ■ (2.12a,b) 

1N m 1 Q > a 


Substitution of equations (2.12) into equations (2.9) yields the following governing 
equations for no impact case 


f n 2 

ah b 


q al I ~~2 

VOX) + « n i 


q = F_ + 1 ; 
» 4 ml m 


(2.13a,b) 


f 2 


h I = <t>p-ta n 


- l 




(2.13c) 
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Figure -2.9. Illustration of different impact regimes. 


Hence, when there is no tooth separation, the system is linear and the mean and 
alternating components of the response, and q^, are uncoupled from each other. 

The transition frequencies from the no impact (linear) regime to the other non-linear 
regimes where tooth separation occurs are found by substituting q a = Id m - *1 i* 110 
the equations (2.13a) and (2.13c) as 



(2.14) 


where f2 tl and ft t2 are the transition frequencies from no impact to single sided 
impact regimes below and above resonance, respectively. 

Case II: Single Sided Impact : Mathematically, the single sided impacts (tooth 
separation without back collision) are observed if 

+ and | t lm _c la| <1, (215) 

As illustrated in Figure 2.9 as case n, the solution remains in the region q(t) > - 1 . 
The describing functions satisfying equation (2.15) are in the form 


N m =1 - 2q 


an_ 

mil 


[Y+ — g( Y— ) ] ; N a =l-i[l + h(Y->] . (2.16a,b) 
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Substitution of equations (2.10a) and (2.10b) into equations (2.16a) and (2.16b) 
gives the describing functions for the single sided impact case: 


N m = i+ ^{M i+ (^) ,i '-> 2 ]-( ir *)} ; <2i7a) 

n mil 

V, 1 2 ( 1_q mn)r, M-llV.. .21 

N » = 2 L* " V 4 Y "> J' <2.l7b) 

Then equations (2.9a), (2.9b) and (2.9c) with the describing functions given by 
equations (2.17a) and (2.17b) define the response of the system in the single sided 
impact region. In this case, it is hard to find closed form expressions for transition 
frequencies from single sided to double sided impact regions. The validity of the 
solutions, obtained by solving equations (2.9) and (2.17) iteratively, should be 
checked. If the solution does not satisfy equation (2.15), then single sided impacts 
are not seen at that particular frequency. 

Case III: Double Sided Impact : The double sided impact case exists if q,,, and 
qg are such that the following conditions are satisfied: 

9a > l 1 “^ml and q a >|l + qj (2.18) 

Figure 2.9 (case IE) shows the double sided impact case at which q a is large enough 
when compared to q^ so that back collision is observed. The describing functions 
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for this case are obtained from equations (2.8a) and (2.8b) under the condition given 
by equations (2. 19): 


N = 1- 


m 


2(tc - 2) 
71(1 am 




1 

^aTTI 



(2.19a,b) 


The solution for the double sided impact case is obtained by solving equations 
(2.9) and (2.19) numerically. The validity of the solution, again, should be examined 
by using the conditions defined by equation (2.18), as it is done in the single sided 
impact case. 

The solutions of all three regimes are combined to obtain the overall frequency 
response of the gear pair. Figure 2.10 illustrates typical q a versus and q,,, versus 

£2 h plots. All three impact regimes are shown on these plots. Also we note the 
frequency region where multiple steady state solutions are obtained. 

2.5. COMPARISON OF TWO SOLUTION METHODS 

First, we validate the approximate analytical solutions of equation (2.5) by 
comparing predictions with the results obtained by digital simulation. Again, two 
example cases (a heavily loaded system and a lightly loaded system) are considered 
and the frequency response curves q^, versus and q,,, versus O h are generated. 

A 

Figure 2.1 1 shows the frequency response for the heavily loaded system with F =2, 
F m =0.1 and £=0.05. Numerical and analytical results agree very well as both predict 

amplitudes and the transition frequencies for case I (no impact) and case II (single 
sided impact) regimes, and both show that the double sided impact solutions do not 
exist. These results demonstrate that the analytical solutions are indeed correct for the 
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frequency 


Figure 2.10. Typical frequency response plots for a gear pair; a) q a versus £2 h , b) q m 
versus 
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heavily loaded system. Now consider the lightly loaded system as the second 
example case with F=0.5, F m =0.1 and £=0.05. From Figure 2.12 we note that 

while case I and case II solutions as yielded by both solution techniques are very 
close to each other, case HI regime is predicted only by the analytical expressions. 
Why is case HI not predicted by the digital simulation? To answer this consider the 
following: 

a) In digital simulation, several sets of q(0) and q(0) must be tried to find all 

of the steady state solutions within a multi-valued region, for instance 

0.6 £ 0. 7 in Figure 2.11 and 0.4 < £0. 7 in Figure 2.12. When q(0) 

h to 

and q(0) are kept constant for each ftj, or if the steady state solution of the previous 
frequency is used as the basis for the initial conditions for the next frequency 
considered, only one of the solutions can be found while missing the other(s). 
Figure 2.13 illustrates this point as the steady state solution of the previous frequency 
is used here for the initial conditions at each Using this procedure, case II 

solutions within the multi-valued regions are missed. Therefore, whenever digital 
simulation is used to solve the equation of motion, dependence of the steady state 
solutions on q(0) and q( 0) must be taken into account to avoid the risk of obtaining 

an incomplete frequency response description. 

b) In the case of Figure 2.12, no initial conditions governed by case III are 
found by the digital simulation within the range of initial conditions - 2 < q(0 ) < 2 
and - 2 < q(0 ) < 2 we have considered. Conversely, the analytical method 
predicts this regime as the issue of initial conditions is irrelevant here. To illustrate 
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this, first we consider the case of Figure 2.3 with F=0.5, F m =0.1, £=0.02 and 
£2 h =0.5. Figure 2.14a shows the domain of attraction governed by case HI within 
the range - 2 < q(0 ) < 2 and - 2 < q(0) < 2 with a q(0) or q(0) increment of 

0.2. Here, almost half of the initial conditions considered define case DI solutions. 
Therefore the chance of having a double sided impact solution at Q h =0.5 is 

considerably high. Now, increase the damping to £=0.03. Once again a case m 
regime is predicted analytically, but the number of initial conditions corresponding to 
case El is not sufficiently high like the case of £=0.02. In the digital simulation, the 
transient solution will converge to case I or case D solution unless the initial condition 
corresponds to one of those shown in Figure 2.14b for case HI. And the limiting 
case is reached when £=0.04, while holding the other parameters the same; now no 
initial conditions corresponding to case HI are found. Obviously there might be 
initial conditions, out of the range considered, which correspond to the analytical 
results. Therefore the existence of case III solutions should be checked numerically 
by searching the entire range of initial conditions as defined by the physical 
considerations of the system. 

2.6. PARAMETRIC STUDIES 

Frequency response amplitudes, transition frequencies and the existence of 
various impact regimes depend on F m , F and £. Therefore, a set of parametric 
studies using analytical solutions of section 2.4 will be presented here; same results 
can be duplicated by the digital simulation technique. First we examine the effect of 

A 

F m and F while holding the damping ratio £ equal to 0.05. Figures 2.15 and 2.16 
compare results for four different F values obtained by varying F ah for a given 
F m =0.1. Both gears maintain complete contact with each other when F is very 
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Figure 2.11. Comparison of theory with numerical simulation; F m =0.1, 
a)q a versus £2 h , b) q m versus £2 h . 
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Figure 2.12. Comparison of theory with numerical simulation; F m =0.1, F =0.5, £=0.05; 
a) versus Q h , b) q^ versus ft h . 







amplitude amplitude 


(a) 



(b) 



Figure 2.13. Comparison of theory with numerical simulation when initial conditions are 
kept constant; F m =0.1, £=0.05; a) F=2, b) F=0.5. 
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Figure 2.14. 
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Domains of attraction for case HI solutions for F m =0.1, F =0.5, £2 h =0.5; 
a)C=0.02, b) ;=0.03. 
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large, say F=10 in Figures 2.15a and 2.16a. Consequently the dynamic system is 
linear and the mean q m and alternating q a components of the torsional motion are 

uncoupled as expected. However, for F =2, no-impact regime (case I) can not cover 
the whole frequency range and a region around the resonant frequency is dictated by 
the single sided impacts (case II). Away from the resonance, over the range given by 
i2|j>0.8 and ^<0.7, solutions are single-valued, whether case I or case II solutions. 
Conversely, near the resonance given by 0.7<£ljj<0.8 dual-valued solutions are seen 

in Figures 2.15b and 2.16b. And sudden discontinuities are observed, i.e. a jump 
up at f2jj=0.8 and a jump down at £2^=0. 7. When F is reduced further, double 

sided impact (case HI) solutions start to appear in Figures 2.15c and 2.16c. Now, the 
frequency range of interest can be divided into the following four regions: i) case I, 
ft h <0.4; ii) case I and case m , 0.4<f2 h <0.6; iii) all three cases, 0.6<£2 h <0.625; and 

A 

iv) only case n, £2 h >0.625. Region (ii) disappears at F >0.5. 

Next, we vary F by changing F m for a given F ah =0.05. This set of 
parameters does not yield any double sided impacts (Case HI). Again, the tooth pair 
does not lose contact when F is large enough, say F^10 and only the linear 

A A 

solutions exist as shown in Figure 2.17a. But when F is lowered to F =4, 2 and 1 , 
the response is non-linear which is composed of cases I and II as shown in Figures 
2.17b-d. 

Next we examine the role of damping ratio £ on the frequency response in Figure 
2.18, given F m =0.1 and F ah =2. Double-sided impacts are found at a low damping 

value £=0.025 as shown in Figure 2.18a. When £ is increased to 0.05, case III 
solutions no longer exist and case I and case II solutions define the frequency 
response completely. The jump-up and jump-down transition frequencies in Figure 
2.18b are distincdy apart. With an increase in £, transition frequencies approach each 
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Figure 2.15. Frequency response q a for different alternating loads; F m =0.1, £=0.05, 
a)F=lO, b) F =2, c) F=l,d) F=0.5. 
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Figure 2.18. Frequency response q a for different damping values; F m 0.1, 
a) £=0.025, b) £=0.05, c) £=0.075, d) £=0.1. 
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other which narrows the dual-valued solution region. At £=0.1, the frequency 
response is single valued and the jump-up and jump-down transitions in Figure 2.18d 
take place almost at the same frequency. Finally, like the linear systems, the damping 
ratio also lowers the amplitudes in the non-linear resonance regime. 

2.7. EXPERIMENTAL VALIDATION 

Our analytical solutions of Section 2.4 will be compared with the experimental 
results of Munro [27] who used a four-square test rig to measure the dynamic 
transmission error x(T) of a spur gear pair. High precision spur gears with 
manufacturing errors much smaller than tooth deflections were selected. Pinion and 
gear were identical with 32 teeth, face width of 12.7 mm and diametral pitch of 4. 
Tooth profile modifications were applied to obtain a minimum (but not zero) e(T) at 
the design load of 3780 N. Other components of the set-up including shafts, 
bearings and casing were made as rigid as possible in order to simulate the 
configuration shown in Figure 2.1. x(T) was measured for a range of gear mesh 
excitation frequencies under different mean loads F m . Some of the key parameters 

were not specified by Munro [27]. For example, it was stated that some additional 
inertias were added to the gears to shift the primary resonant frequency within the 
operational speed range, but the specific values of such inertias were not given. It 
was also reported that the damping ratio £ varies with load "in a random manner". 

Also, backlash was not measured or reported. Therefore, in our study, we estimate 
the damping ratio (£=0.017.5) and the resonant frequency co n by considering the 

design load case at which only linear behavior is seen; time-invariant mesh stiffness is 
assumed in the model. The same value of £ is used at each discrete load F m and a 

backlash value of 2b=0. 1 mm is assumed in our model. 
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Figures 2.19 to 2.22 compare measured and predicted frequency response of 
the dynamic transmission error x(t ) , on a peak to peak basis, at the design load and 
at 1/4, 1/2 and 3/4 of the design load, respectively. As shown in Figure 2.19, both 
analytical and experimental results indicate that teeth maintain contact at the design 
load and hence the system behaves as a linear system in spite of the backlash. This is 
because the static transmission error e(F) is minimum at the design load which 
results in a large force ratio F m /F a , say 50. However, when the mean load F m is 

lowered to 3/4 of the design load which corresponds to a larger static transmission 
error e(t), tooth separation takes place. Consequently we note a jump in the 
frequency response as shown in Figure 2.20 for both analysis and experiment. This 
jump phenomenon is more noticeable at 1/2 and 1/4 of the design loads as shown in 
Figures 2.21 and 2.22, respectively. Our theory matches very well with Munro's 
experimental results [27] in spite of the lack of knowledge of some relevant 
parameters. 

As the second example case, the experimental results of Kubo [28] are 
considered and compared with our theory. Experimental results and the relevant 
system parameters are extracted from a recent paper by Ozguven and Houser [24] and 
the excitation e(D is calculated using a spur gear elastic model [21]. Kubo 
designed a four-square spur gear test rig which was heavily damped (C“0.1). He 
measured dynamic root stresses cr d and then estimated the dynamic factor D s as 
D c =Oa/o where O c is the static tooth root stress. However, in several other studies 
[16] the dynamic factor is defined as the ratio of the dynamic mesh load F d to static 
mesh load F m , given by D L = F A / F m = F F m = (2£q(t) + f(q(t) ) / F m . 
Note that D L is equal to D s when the change of the moment arm due to a change in 
contact point is neglected. Figure 2.23 shows the envelope of dynamic factor D s 
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gure 2.19. Comparison of the theory with Munro's [27] experiments at design loau. 
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Figure 2.21 . Comparison of theory with Munro's 
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Comparison of theory with experimental results of Kubo (extracted from 
reference [24]). 
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measured for eight different teeth pairs and spectrum computed using the 
analytical solutions of Section 2.4. Here, at most of the frequencies, our predictions 
are not within the measured envelope, but the transition frequency and the amplitude 
at the jump discontinuity are predicted accurately. There might be several reasons for 
the discrepancy including the usage of instead of D<j, the validity of the 

computer model used to calculate e(t ) , insufficient knowledge of some system 
parameters such as profile modifications, and the assumptions made in developing 
our theory such as the time invariant mesh stiffness. Therefore the experimental data 
is not exactly analogous to our analytical predictions even though a satisfactory 
agreement is found. Obviously our theory should be refined in order to obtain a 
better agreement with Kubo's data; for instance, time-varying mesh stiffiiess may be 
included. 

2.8. COMPARISON OF EXCITATIONS 
2.8.1. Internal versus External Excitation 

First, we assume that only one type of excitation exists at a time and compare 
the frequency response characteristics of the system due to the internal static 
transmission error sinusiodal excitation Fj(t) given by equation (2.5) with the external 
sinusoidal excitation F e (t) given by equation (2.4). The comparison is based on the 
analytical solutions which are constructed in Section 2.4 for equation (2.5) and in 
Reference [33] for equation (2.4). In the case of internal excitation, the amplitude of 

the alternating force has a term which makes the alternating force amplitude 

2 

frequency dependent. This amplitude F h is smaller than F ah for f2 h <l and 
greater than F^ when fl h >l. Hence the overall alternating force amplitude ratio 
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^ah^h var ‘ es even thou gh F = F m / F^ is kept constant. In the 

case of external excitation, the force amplitude ratio F=F m /F aT is frequency 
independent. 

Figures 2.24 and 2.25 compare the frequency responses for four values of F 
given F m =0.1, f2 h =Qj=n and ^=0.05. When F is sufficiently large, say F £10, 

only case I solutions corresponding to the no-impact (case I) exist. Accordingly both 
excitations result in the linear system responses which are close to each other for 
ft<l; but for £2>1 there is a difference which grows with increasing frequency as 
shown in Figure 2.24a. When F is reduced to 2, both case I and case II regimes 
exist; the difference between two excitations is again significant at higher Q as shown 
in Figure 2.24b. The transition frequencies also differ, and a larger range of dual- 
valued solutions is seen for equation (2.5). Case HI solutions are witnessed at lower 

A 

values of F (1 and 0.5) as shown in Figures 2.25a and 2.25b. While equation (2.5) 
always has a case I regime at low Cl, equation (2.4) does not produce case I at 

A 

F =0.5. Another important observation from Figures 2.24 and 2.25 is that up to two 
steady state solutions are seen in the external torque excitation case; conversely as 
many as three solutions are found for the geared system excited by the static 
transmission error. 

2.8.2. Periodic or Combined Excitation 

The approximate analytical solutions given in Section 2.4 are constructed only 
for a single harmonic internal excitation term e(F) or Fj(t). However, in the real 
geared systems, e(t ) or Fj(t) is a periodic which can be represented by a fourier 
series of the fundamental frequency £2 h . Therefore, it is necessary to consider the 
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relative displacement, x relative displacement, 


(a) 



(b) 



Figure 2.24. Comparison of the frequency response due to the external torque excitation 
given by Comparin and Singh [33] and the static transmission error excitation; 
£=0.05, F m =0.1, n h =Q x and a) F=10, b) F=2. 
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frequency 



frequency 


Figure 2.25. Comparison of the frequency response due to the external torque excitation 
given by Comparin and Singh [33] and the static transmission error excitation; 
£=0.05, F m =0.1, Q h =Q x and a) F=l, b) F=0.5. 
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higher harmonics of Fj(t) besides the fundamental component which is already 
included in the analysis. On the other hand, both external F e (t) and internal Fj(t) 
excitations may exist simultaneously. These two cases require that F(t) in equation 
(2.3) be reformulated as follows 

k 

F(t) = F m + F e (t) + F.(t)= F m + £ F aT .sin(ja T t+4> Tj ) 

j= 1 

+ I (j« b ) 2 F ahj sin < h' + +hj> (2.20) 

j = 1 

In order to construct analytical solutions of this equation, we must investigate 
the applicability of the principle of superposition, which would consider first each 
excitation separately and then superimpose the corresponding responses to generate 
the overall frequency response. Two cases are considered and the analytical results 
based on the principle of superposition are verified through digital simulation. 

Consider only the static transmission error excitation Fj(t) in equation (2.20) 

3 2 

with 3 harmonics, i.e. F e (t) and F(t) = F m + X (j^h) Fahj + < l , hj)- The 

j=l 

amplitudes of the harmonics are selected to be F^^O.05, F ah? =0.02 and F ah ^=0.01, 
and the mean force F m =0.1. First each excitation component is considered separately 
without paying any attention to the phase angles and the frequency responses are 
obtained using digital simulation as shown in Figure 2.26. Second, all there 
excitations are included simultaneously and the overall response is compared with the 
previous solutions in Figure 2.26. The frequency response due to each harmonic of 
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the excitation is very close to the overall solution around the resonance peaks 
governed by that harmonic. When all of the solutions corresponding to each 
harmonic are added algebraically by also considering the phase angles, the overall 
response improves as shown in Figure 2.27. This response is in the form 


qa(0 - X q a jsin(ft h t + <f>hj) where q a j is the alternating displacement when only 


the j-th harmonic is considered. Figure 2.27 suggests that the principle of 
superposition can be applied to a gear pair with backlash provided that the excitation 

frequencies are sufficiently far apart from each other. In the case of periodic static 
transmission excitation, this is valid as all excitation frequencies are at least apart. 

Now, the principle of superposition, which has been already verified by digital 
simulation, can be used to obtain the approximate analytical solutions per Section 2.4 
when the periodic forcing function is considered. Figure 2.28 shows the analytical 
frequency response curve for F m =0.1, F ahl =0.05, F ah2 =0.02, F ah3 =0.01 and 
£=0.05. Here the jump discontinuity is seen only at the peak Q hl governed by F^, 
since F m /F ah2 and F m /F ah3 are both sufficiently high so that no tooth separation 
occurs at and fi^ 3 , respectively. However, the jump phenomenon can also be 
seen at the higher harmonics depending on the force ratios ^/F^j, j* 1 and £ in 
accordance with the results of Section 2.6. The same concept can be applied to the 
superharmonic components of the external torque excitation of equation (2.20) 
provided Fahj 0. 

The principle of superposition is now extended to the case where both internal 
Fj(t) and external F e (t) sinusoidal excitations exist simultaneously as given by 
equation (2.3) provided that the excitation frequencies f2 h and Qy are not close to 
each other. In a real geared system, fijj is much higher than fly which implies that 
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Figure 2.26. Comparison of the frequency response q a due to first there harmonics of the 

static trancmission error for the cases when the harmonics are considered 
simultaneously and separately; F m =0.1, £=0.05, F^^O.05, F ah2 =0.02 and 
F^^.Ol. These results are obtained using digital simulation. 
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first 3 harmonics simultaneously 
principle of superposition 



frequency 


Figure 2.27. Comparison of digital simulation results obtained by using the principle of 
superposition and by considering all three harmonic excitations applied 
simultaneously; F m =0.1, £=0.05, F^^O.05, F ah2 =0.02 and F^sO.Ol. In 

this case, phase angles are included in the superposition. 
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amplitude 



Figure 2.28. Comparison of the frequency response q a due to first three harmonics of the 
static transmission error. This curve is based on the principle of 
superposition applied to the analytical solutions which have been constructed 
by the method of harmonic balance; F m =0.1, £=0.05, F^^.05, F ah? =0.02 
and F aM =0.01. 
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frequency 


Figure 2.29. Frequency response q a due to both the external torque and the static 
transmission error excitations. This curve is based on the principle of 
superposition applied to the analytical solutions which have been constructed 
by the method of harmonic balance; F m =0.1, £=0.05, F^=0.05, F ^0.05, 

n h =2n T . 
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the principle of superposition should be suitable for this case. Figure 2.29 illustrates 
the frequency response solutions due to the fundamental harmonic component of both 
excitations for F m =0.1, F ahl =F tT1 =0.05, C=0.05 and by assuming that fl h =2Q T . 
Here, the jump discontinuity is seen at two frequencies. However, when £2 h = H T 

or when one of the superharmonic peaks of Fj(t) coincides with the resonant peak 
governed by F e (t), the principle of superposition will no longer be applicable. In 
such cases, we will use the digital simulation as the analytical interpretation of these 
cases is yet to be explored. 

2.9, CONCLUSION 

This analytical study on the non-linear dynamics of a spur gear pair with 
backlash as excited by the static transmission error has made a number of 
contributions to the state of the art. First, difficulties associated with the digital 
simulation technique have been resolved as multiple steady state solutions at a given 
frequency can be found provided the entire initial condition map is searched. Second, 
new frequency response solutions for the gear pair have been constructed using the 
method of harmonic balance. Third, our mathematical models have been validated as 
these compare well with two previous experimental studies and the key parameters 
such as the mean load, mean to alternating force ratio, damping and backlash have 
been identified. Forth, the chaotic and subharmonic resonances are observed if the 
mean load is too small for a lightly damped system. Fifth, mathematical conditions 
for tooth separation and back collision have been established which are compatible 
with available measured data. Sixth, the periodic transmission error excitation case 
has been analyzed using the method of harmonic balance in conjunction with the 
principle of superposition. Finally, on a more fundamental note, our study enriches 
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the current literature on the clearance non-linearity or vibro-impact systems as the 
governing non-linear differential equation is different from the conventional single 
degree of freedom system formulation. 



CHAPTER III 


NON-LINEAR DYNAMIC ANALYSIS A GEARED ROTOR- 
BEARING SYSTEM WITH MULTIPLE CLEARANCES 

3.1. INTRODUCTION 

Mathematical modeling of geared rotor-bearing systems, being an essential step in 
designing quiet and reliable power transmissions, has been the subject of numerous 
studies over the past few decades. Most proposed dynamic models, as reviewed by 
Ozguven and Houser [14], are essentially linear. However in several cases, it has been 
experimentally shown that the geared systems exhibit non-linear behavior [27-29,35]. 
For instance, vibro-impacts are observed in a lightly loaded transmission with gear 
backlash or loose bearings [30,31,34], Accordingly, we must develop non-linear 
mathematical model of the geared system - this is main focus of this chapter with 
emphasis on the clearance type non-linearity in gears and rolling element bearings. In 
Chapter II we had considered the single degree of freedom non-linear model of a spur 
gear pair with backlash and investigated the effects of system parameters on the 
vibrations and chaos excited by the static transmission error. Applicability of the 
harmonic balance method and digital simulation technique to the solution of the steady 
state response has been demonstrated, and difficulties associated with the digital 
simulation technique when applied to such systems governed by stiff non-linear 
differential equations have also been resolved in Chapter II. 

Although there is a vast body of literature which considers a single degree of 
freedom system with clearances, as reviewed earlier in Chapter D, studies on the multi 
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degree of freedom vibro-impact systems are very limited. For example, Galhoud et.al. 

[52] considered a two degree of freedom translational system with a gap and found the 
forced harmonic response using the piecewise linear technique. Winter and Kojima 

[53] also used the same technique to study the geared systems with backlash. 
However, it should be noted that the piecewise linear technique can not predict several 
non-linear phenomena such as subharmonic and chaotic responses since it employs the 
assumption that both impact and no-impact regimes "are repeated in an identical manner 
once every period of excitation" [52]. Kucukay [35] developed an eight degree of 
freedom model of a helical gear pair with backlash to include the rocking and axial 
motions of rigid shafts and radial deflections of linear bearings. Lin et.al. [38] included 
motor and load inertias in a three degree of freedom torsional model. Both of these 
studies employed the digital simulation technique, but did not consider a number of 
issues of primary concern in non-linear system such as the existence of multiple steady 
state solutions, their dependence on initial conditions, subharmonic, quasi-periodic and 
chaotic responses, etc. A few of these issues have been addressed by Singh et.al [34] 
on the gear rattle problem. Comparin and Singh [54] have also used the digital 
simulation, analog simulation and harmonic balance method to analyze coupled impact 
pairs assuming that the modes are "weakly coupled", which allows the system to be 
represented by a combination of single impact pairs. They included only the low 
frequency external torque excitation and found the steady state frequency response at 
the primary resonance. This solution has then been used to analyze the neutral gear 
problem in more detail [55]. 

The other two groups of studies which consider multi-degree of freedom systems 
with continuous non-linearities [56-61] and periodic excitations due to mesh stiffness 
variations [29,62], will not be addressed here since their responses are significantly 
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different from the clearance non-linearities as discussed by Comparin and Singh [33]. 
The literature on chaotic vibrations and rolling element bearings will be discussed later. 

3.2. PROBLEM FORMULATION 
3.2.1. Scope 

Figure 3.1 shows the generic geared system considered in this study. It consists 
of a spur gear pair mounted on flexible shafts which are supported by rolling element 
bearings assembled in a rigid gear box. Since the shafts and bearings are compliant, 
our single degree of freedom model of spur gear pair which assumes fixed gear centers 
is obviously not suitable. Instead, a three degree of freedom non-linear model as 
shown in Figure 3.2a is considered. It includes equivalent stiffness and damping 
elements representing the shaft and the bearings. The corresponding linear model has 
been found to be sufficiently accurate when compared with a finite element model for 
eigen solutions provided that the gear dynamic response (mesh force, dynamic 
transmission error, etc.) is of major concern. Gear backlash and radial clearances in 
bearings are defined analytically. An approximate non-linear bearing model is also 
proposed. Applicability of both analytical and numerical solution techniques to this 
problem is investigated. Several key issues such as non-linear modal interactions and 
differences between internal static transmission error excitation and external torque 
excitation are discussed. Parametric studies are conducted to understand the effect of 
system parameters such as bearing stiffness k^j to mesh stiffness kjj ratio 

k = k bi / k h , alternating to mean force ratio Fj/Fjj, and radial bearing preload to mean 
force ratio Fb/Fm on the non-linear frequency response. A criterion is used to classify 
the steady state solutions, and the conditions for chaotic, quasi-periodic and 
subharmonic steady state solutions are determined. Two typical routes to chaos 
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rolling 

element 

bearing; 


Figure 3.1 . Generic geared rotor-bearing system. 



observed on a geared system are also identified. Finally, our formulation is verified by 
comparing predictions with measurements [27]. 

3.2.2. Physical Model and Assumptions 

The three degree of freedom non-linear model of the geared rotor system with 
gear inertias Igj and Ig 2 , gear masses mgj and nig 2 , and base circle diameters dgi and 
dg 2 , as shown in Figure 3.2a, is considered here. The gear mesh is described by a 
non-linear displacement function f b and viscous damping c b . Friction forces at the 
mesh point are assumed to be negligible [18]. Thus the transverse vibrations in the 
pressure line direction are uncoupled from the vibrations in the direction perpendicular 
to the pressure line. Bearings and the shafts that support the gears are modeled by 
equivalent elements with viscous damping coefficients c bl and c b 2 and non-linear 
springs defined by force-displacement functions f b j and f b 2- The effects of the prune 
mover and the load inertias are not considered assuming that these inertial elements are 
connected to the gear box through soft torsional couplings. Further, we assume that 
the system is symmetric about the plane of the gears and the axial motion (parallel to the 
shafts) is negligible. Like the spur gear pair model of Chapter II, both low frequency 
external excitation due to torque fluctuations and high frequency internal excitations due 
to the static transmission error e(t) are considered in the formulation. Input torque 

fluctuation is included, but the output torque is assumed to be constant, i.e. 
Tgi(t) = T glm +Tgj a (t) and T g2 (t) = T g2m . External radial preloads ^ and i\, 2 are 

also applied to both rolling element bearings. 
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3.2.3. Equations of Motion 

Equations of coupled transverse-torsional motion of the non-linear geared rotor- 
bearing system shown in Figure 3.2a are 

m gi ygi + c bi ygi + c h (*' + ygi - yg2 - «') + * b i f bi(y g i ) 

+k h f h (x + y g i-y g 2 -e) = -^,i (3.1a) 

m g2 yg2 +c b2 yg2 _c b(^ , + ygl ~yg2 _ ® / ) + k b2 f b2(yg2) 

-k h f h (x + y gl -y g 2 -e) = ^ ) 2 (3.1b) 

m cl *" + c h (r + ygi - y g 2 - e') + k h f h (x + y gl - y g2 - e) = F m + F aT (t) (3.1c) 


s(i) = ^e s ,(i)-iie g2 (i) ; 


m cl — / 2 

' d 2 


(3.1d,e) 


gl u g2 


41 


gl 4I g2 


l m 


^^glm _ ^^g2m p /-\ _ m cl^gla(0 
d g i d g2 ; aTU_ 2I gl 


(3.1f,g) 


Here, ( )' means derivative with respect to time t , y gi and 0gj are the transverse and 

torsional displacements of the i-th gear (i=l,2), rr^j is the equivalent gear pair mass, 
F m is the average force transmitted through the gear mesh and F aT (t) is the fluctuating 

force related to the external input torque excitation. Equations (3.1a), (3.1b) and (3.1c) 
are simplified further by defining below a new variable p(t) which is the difference 
between the dynamic transmission error and the static transmission error e(t). 
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(b) 


Figure 3.2. Non-linear models of the geared rotor-bearing system; a) three degree of 
freedom model, b) two degree of freedom model. 
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p(t) = x(t) + y gl (f) - y g2 (t) - e(t) ; 


(3.2a) 


m gi ® ® ygi(t) cj,j o cj, fygiO)] 

0 m g2 0 yg2(t) + 0 C b2 -c h y' 2 (t)l 

-m cl m cl m cl J[ p^(t) J L° 0 c h J p'(t)J 

k bl 0 k h 1 f bl<ygl) -^,1 

+ 0 k b2 -k h - f b2 (y g2 )> = < fy 2 \ (3.2b) 

-° 0 k hJ[fh(P) J [pm “m cl e"(t)+F aT (t> 


A dimensionless form of equation (3.2b) is obtained by letting 
y g i(0 = ygi(t)/b c , p(t) = p(t)/b c , ® n = V k h/m c i, ©w = ^/k bi / mgi (i=l,2) and 
t = © n t where b c is the characteristic length. Here, we consider harmonic excitation for 
both e(t) and F aT (t) as e(t) = esin(ft b t + <J> b ), F aT (t) = ^aT sin(ft r t + <(>t) where fl b 

and are the fundamental excitation frequencies of internal displacement and external 
torque fluctuations, respectively . Further, define dimensionless excitation frequencies 
= / ©n ©td = / ©n to yield the following dimensionless governing 

equation of motion. 

'i o ol|y g i(t)l fa, o Ci 3 lfy g i(t) 

0 1 0 * y g2 (t) *+2 0 C22 -C23 *yg2( t )‘ 

. - -1 1 . p(t) , 0 C33. . p(t) 

Kh 0 K 13 f b i(y g i) 

+ 0 K 22 -k 23 < f b 2 (y g2 ) ={F(t)} ; (3.3a) 

-° 0 1 J. f h(P) J 
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sin(£2 T t + <t>T>;(3.3b) 


{F(t)} = {F} m + {F(t)}j + {F(t)} e 


-hi 


0 


0 

hi 

^ + < 

0 

>sin(£2 h t + < hi) + ’ 

0 - 

, F m , 


Fah^b 


F aT. 


Cii = 


— r,, = — E b_ 


’ ^3 “ ■ I1 r ” - i=1 - 2 i C33 = 

igjiu a 2mgjto n 


_ c h . 


(3.3c, d,e) 


K i3 = — ; Fi»= — &i-y. i=1.2; 

(D n m gi m gi^c w n 


(3.3f,g,h) 


1 m 

m ~ 


m cl^c®n 


r e F _ F aT 
2> F ah = — ; F aT - _ “ 


b c ’ 


m cl b c ®n 


(3.3i,jjc) 


where Fbi (i=l,2) and F m , are the dimensionless components of the mean force vector 
{F} m , and F aT and F^ pertain to the alternating external excitation {F(t)} e and internal 
excitation {F(t))j force vectors, respectively. 

3.2.4. Modeling of Non-linearities 

The non-linear displacement functions fbi(ygi), i=l, 2 and f^Cp) in equation 
(3.3a), which represent the bearing radial and gear mesh stiffnesses, respectively, 
should be defined explicitly before solving the non-linear equations. Here, fb(p) is, 
defined as a clearance type dead space function with a backlash 2bh and linear time- 
invariant mesh stiffness of one, in the dimensionless form 
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I p — T p 

I b.. b. 


,(p) - v' • 




h - P " ^ 


(3.4) 


Or. 


P < 


For the i-th roll in;: element bearing, radial three F v i versus displacement Vg. 
relationship unde- the static loading condition is defined x< j 03-65 i 


\H 


l K «' ^ \>Vi COS(t i 


■ J cos a, ; 


egi " " Di 


My gi 1 


F bi <y» ^ F bi 


H „ -r 

i. 1 !|ye«K-«w;o- J - bj n j 


< -~b 


hi 


(3.5) 


where kq h; the nnei contact stiffness, ft T is the angular position of the r-th rolling 
element in contact, 2b- 0l is the radial clearance of the i-th hearing., n is the power of the 
non-linear force displacement relationship (n=1.5 for hall hearings and n~10/9 for 
roller bearings) and H is the total number of rolling elements in contact under loaded 
conditions. Now the dimensionless bearing displacement function equation 

(3.3a) is obtained from equation (3.5) as 
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Hf 

I 

I T=1 


Dk; 

y^cosaj — — cosa r ; 
b c J 


y . > 

yg « > b 
U C 


f bi(ygi)- y b n -{ 

lt u u c 


0 ; 


H (> 

-I 

r=l 


_ b bL <v . <^bi 

k y «‘ k 


(3.6) 


y gi 


cosa r - 


b bi 




cosa r ; 


y g< 


<- 


’bi 


J c J 


Figure 3.3 shows the non-linear function fgj(y g j) for a roller bearing with n=10/9, 
k t =lxl0^ N/mlO/9, total number of rollers Z=15 and bjj=0.01 mm. In Figure 3.3, we 
note almost a linear r n1< »tionship for large displacements, say y g i>-3b] ;> j/b c . Since the 
degree of non-linearity is not very significant, equation (3.6) can be approximated by a 
piece-wise linear function, similar to fj,(p) given by equation (3.4), in order to simplify 
the analysis considerably. Figure 3.3 also illustrates two linear approximations A and 
B beyond the clearance for fbi(ygj) in the form 

f bi(y g i) = 


Note that both approximations A and B differ in clearances b^ and b^jg but have the 
same slope as the exact bearing stiffness curve for ygi>3bjjj/b c , i=l,2, which is unity in 
the dimensionless form. The validity of these approximations will be given later. 


y g < 

0 ; 

y gi 


b bi . 
b r ’ 


b bi . 


b bi 

** > i L 

U C 

hfeL < v < --bi 
b c <y « i< b e - 


(3.7) 


y gi<- 


>bi 
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3.3. CORRESPONDING LINEAR MODEL 

As a limiting case, equations of motion of the corresponding linear system are 
obtained by substituting fbi(ygi ) = ygi , for i=l and 2, and fh(p) = p in equation (3.3a) 


1 0 o' 

y g i(0' 


Cl 1 0 Cl 3 

y g i(t) 


0 1 0 

y g2 (t) 

> + 2 

0 C 22 - C 23 

y g2 (0 

► 

-1 1 1 

,p(t) . 


Lo 

0 C 33 J 

p(t) 

- 




’*11 0 

* 13 ] 

ygiD) 



+ 

0 *22 

— *23 ' 

y g2 (t) ► 




0 0 

*33. 

. P(0 . 


= {F(.)} 


(3.8a) 


or in the matrix form with ( q(t) } as the displacement vector 

[M]{q(t)| +[C]{q(t)} +[K]{q(t)| = (F(t)} (3.8b) 

where the mass [M], damping [C] and stiffness [K] matrices are all positive definite. 
These matrices are asymmetric due to y g j and y g 2 terms in the last row of matrix 

equation (3.8a). For this linear system, the mean and alternating components of the 
motion can be separated by letting y g i(t)=y g i m +y g i a (t), i=1 »2 and p(t)=p m +p a (t). 
Hence, equation (3.8a) is rearranged in terms of the alternating motion as 

[M]{q a (t» +[C](q a (t)} +[K]{q a (t)} = {F(t))i +{F(t)} e (3.9) 

The natural frequencies co r and the modes {y r } are calculated by considering the 
corresponding eigen-value problem. The forced harmonic vibration response is then 
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obtained by the modal expansion technique in the following form 


3 T 

(q a (0}=X r — T - J - 2 Vr) 1 F ah ngsin(n h t + «> h ) 

r=l [<®r -^h) + j^h c rj 

+ Z t (VrllVr) » F a -pt2j sin(ft T t + <j> j ); j = V— I (3.10a) 

i=l |^io r -£Tj') + j^T c r 


c r =r^-(V r > T [C]( V s } ; 5rs = { 1; " S (3.10b) 

°rs [0; r*s 


Table 3.1 shows the natural frequencies of the three degree of freedom linear 
system given by equation (3.9) and the ones found by finite element method for three 
different k = k bi /k h values with m gi =l kg, I gi =0.0008 kg/m 2 , d gi =0.08 m, 

Cji=0.01, Cj3=0.0125, ^ 33 = 0 . 05, i=l, 2 and kh=2xl0® N/m. As evident from Table 
3.1, the three degree of freedom and finite element models [15] result in virtually the 
same natural frequencies. This concludes that the three degree of freedom model is 
indeed suitable for the geared-rotor bearing system. The first \|/j and third \j/ m natural 
modes are coupled transverse-torsional modes while is purely transverse type [15]. 
The second natural mode y n is not excited by e(t) in this particular case, since the gear 
ratio y g =d g 2 /d g ] is one for a symmetric pair. Therefore only two peaks should exist in 
the frequency response spectrum which will be presented in Section 3.5. 
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Table 3.1 


Natural frequencies of the corresponding three degree of freedom (3- 
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3.4. TWO DECREE OF FREEDOM SYSTEM STUDIES 
3.4.1. Equations of Motion 

As a first example, we reduce the three degree of freedom transverse-torsional 
semi-definite model to a two degree of freedom non-linear model. This is obtained by 
clamping one of the gear centers as shown in Figure 3.2b, i.e. one of the transverse 
displacements is assumed to be zero, say ygj=0 which physically corresponds to a 
system with one gear gl mounted on a rigid shaft which is supported by very stiff 
bearings while the other gear g2 is assembled with compliant lumped shaft and 
bearings. Equations of motion of the reduced order system are obtained from equation 
(3.3a) by letting y g i=0. 


1 

1 


f 82(,, Ui 

II p«> J 


C22 

0 





f b2(y g 2>l 
*h(p) J 


= {F(t)}. 


(3.11) 


The intent here is to simplify the physical system in order to investigate several 
key issues in depth. Specifically, the objectives are to: a) justify using approximate 
bearing models proposed in Section 3.2.4 by comparing the frequency responses 
excited by c(t), b) show applicability of the harmonic balance method (HBM) to solve 
non-linear system equations and compare its predictions with the results yielded by the 
digital simulation technique, c) study interactions between system non-linearities for 
both weakly and strongly coupled modes, and d) compare the steady state frequency 
response spectra due to internal static transmission error excitation { F(t) )j and external 
torque excitation (F(t)} e . 
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3.4.2. Solution 

Both analytical and numerical solution techniques which we have used 
successfully for the spur gear pair problem are employed again to solve equation 
(3.11). First, approximate analytical solutions are constructed using the harmonic 
balance method [33,51] with only the bearing non-linearity assuming no gear backlash, 
i.e. bh=0 or fh(p)=P- For the harmonic excitation given by equation (3.3) , we assume 

P(t) = Pm +PaSin(fth t + < l>p) I y g 2( l > = y g 2m + y g 2a sin ( ti h t + < t > g2> (3.12a,b) 


where subscripts m and a represent mean and alternating components of the steady state 
response, and <|>p and <|)g 2 are the phase angles. The non-linear bearing function is 

approximated as 


f b 2 (y g 2 > = N m y g 2 ni + N a y g 2 a sin(Q h t + <t> g2 ). (3.13) 

where N m =N m (y g2m ,y g 2a) and N a =N a ( yg2m ,y g 2a) are the describing functions which 
are given in Chapter II; these need to be defined for each unpact regime (no impact, 
single sided impact, double sided impact). Governing frequency response equations 
are obtained by substituting equations (3.12) and (3.13) into equation (3.11) and 
equating the coefficients of like harmonics: 


[A]+(2C 22 ^h) 2 ]^hFah 

Va : 


(3.14a) 
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A = (-flh K 23 A 1 + 4 ^h^22C23 + A 2 A ^ + 4ii hC22 A 2 ^ 


+(— 2Qh^22 K 23 _ 2ti hC23^1 + 2 ^hC33-^ +8 ^hC22C33) 2 i (3.14b) 


Aj — K 22 N a ^b’ -^2 1 » Pm — ^m ’ 


(3.14c,d,e) 


y g 2a = 


*23 + ( 2 ^hC 23) 2 

( k 23 A 1 “ 4 ^b^22C23) 2 + ^ 2ii bC22 K 23 + 2ft hC23 A l 7 



(3.140 


.. _ Fb 2 + * 23 Pm 

yg2m " K 22 N m 


(3.1 4g) 


The overall frequency response is obtained by solving equation (3.14) numerically for 
each impact regime separately. 

Second, the digital simulation technique is used to solve equation (3.1 1) for no 
gear backlash. A 5th-6th order variable step size Runge Kutta numerical integration 
algorithm [50], which has already been employed successfully in solving similar 
problems, is used here. Figure 3.4 compares frequency response spectra obtained 
using the harmonic balance method and digital simulation. Although both methods are 
in a very close agreement for the case considered in Figure 3.4, one should be aware of 
the following: (i) several problems may exist in the application of the numerical 
integration method to the clearance type problems, hence caution must be exercised, 
and (ii) the harmonic balance method is incapable of predicting chaotic and 
subharmonic responses. Accordingly for further analysis, we will use only the digital 
simulation technique. 
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Figure 3.4. Comparison of harmonic balance method and digital simulation technique; 

F m =1.0, F ah =0.5, F b2 =0, b h =(), b b2 =b c , ^3=0.05, C 2 3=0 0125, 
^ 22 = 0 . 01 , k = 1 ; a) y g2a versus Q h , b) p a versus 
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3.4.3. Validation of Non-linear Bearing Model 

We examine the validity of two approximations for fbi(ygi), as suggested in 
Section 3.2.4 and shown in Figure 3.3, by assuming bj,=0 and F a j=0. Equation 
(3.11) is solved using the digital simulation technique, given F m =1.0, F a j 1 =0.05 , 

f r b2 = ®» C22 ss ® 01, £33=0.05, £23=0 0125 and K22=K23=0.25, for three bearing models 
(exact, approximations A: bb2 A =bb2=b c , approximation B: bj,2 B =1.63b|,2=1.63b c ). 

Figures 3.5a and 3.5b compare yg2 a versus Q|, and p a versus spectra for all three 
bearing models. Both approximations are in a very close agreement with the exact 
bearing model over the entire frequency range. The frequency and the amplitude of the 
jump at the first peak is predicted accurately by both approximate models. Hence we 
conclude that the approximate piecewise linear model can be used for a bearing without 
losing any accuracy and the dynamic behavior is not sensitive to the minor variations in 
bearing clearances b^. 

3.4.4. Non-linear Modal Interactions 

Although both natural modes of the corresponding linear model are transverse- 
torsional coupled modes, first mode \|fj is more dependent on transverse vibrations than 

the second peak. Accordingly, it is expected that bearing non-linearities should affect 
the first peak directly while gear backlash should dominate the second peak in Figure 

3.4. The validity of this claim depends on the nature of non-linearities and system 
parameter values, as illustrated below. 

First, a set of parameters is selected so that natural frequencies coj and (Ojj of the 
corresponding linear system are far away from each other: coj =0.437 and Ojp 1.144 
for k 22 = 0 . 25 , K 23 = 0 . 25 , ^22=0. 01, ^ 23=0 0125 , and £33=0. 05 . Frequency response 
characteristics of the system with a forcing function, F m =1.0 and F^O.5 are evaluated 
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Figure 3.5. Comparison of frequency responses obtained by using the exact and 
approximate bearing displacement functions given in Figure 3.3, 
F m =l .0, F ah =0.5, F b 2=0, b h =0, b b2 =b c , ; 33 =0.05, £ 23 =0.0125, 
; 22 =0.01, k = 1; a) y g2a versus Q h , b) p a versus Q h . 
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for the following three cases: a) no gear backlash, b h =0, non-linear bearing, b b 2 “b c ; b) 
gear backlash, bjj=10b c , linear bearing, b^ 2 = 0 i ami c) gear backlash and non-linear 
bearing with b h =10b b 2=10b c . In the last tw F o cases, gear backlash is taken to be 10 
times the bearing clearance as it is a typical order of magnitude value for many geared 
rotor systems supported by roller bearings. Figures 3.6 and 3-7 show the frequency 
responses, y g2a versus £2 b and p a versus respectively, for all three cases given 
above; also superimposed is the linear frequency response curve for comparison 
purposes. As shown in Figure 3.6a and 3.7a, a jump discontinuity is observed at the 
first peak only for case a. In the case with only gear backlash non-linearity, jump is 
seen at the second peak while the first peak is continuous as shown in Figures 3.6b and 
3.7b; this supports the claim of torsional mode dominance of the second peak. Finally, 
when both non-linearities are included simultaneously, both peaks become non-linear 
as shown in Figures 3.6c and 3.7c. For such a system, natural modes are essentially 
’uncoupled" and the interactions between component non-linearities are negligible In 
summary, for case a of Figures 3.6a and 3.7a, a linear model can be used to predict the 
frequency response beyond (0j , say £2 b > 0.6. Similarly, a linear model could be 
suitable for case b for ii b <0.8 as shown in Figures 3 6b and 3.7c But for case c, a 
linear model is suitable only far away from the resonances 

Second, choose a dataset, say k 22 =0.64, k 2 3=0.10, C 22 =0.01, C 23 =0.0125, and 
£ 33 = 0 . 05 , such that two natural frequencies for the corresponding linear system are 
brought closer to each other: copO.720, 0^=1. 101. We again consider the same cases 
a, b and c and compare the frequency response for F m =2.0 and F^O 5 in Figures 3.8 
and 3.9. Here, each non-linearity affects both modes, resulting in jumps at each peak. 
Therefore, for this set of parameters, modes are considered to be "coupled" because of 
dynamic interactions. Accordingly both non-linearities must be included 
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Figure 3.6. y g 2 a versus plots for the case of uncoupled modes; a) no backlash, 
non-linear bearings, b) backlash, linear bearings, c) backlash, non-linear 
bearings; corresponding linear response. 
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a) 


5.0 



frequency 


Figure 3.7. p a versus C2|j plots for the case of uncoupled modes; a) no backlash, non- 
linear bearings, b) backlash, linear bearings, c) backlash, non-linear 
bearings; corresponding linear response. 
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Figure 3.8. yg2 a versus plots for the case of coupled modes; a) no backlash, non- 
linear bearings, b) backlash, linear bearings, c) backlash, non-linear 
bearings; corresponding linear response. 
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Figure 3.9. p a versus ^2^ plots for the case of coupled modes; a) no backlash, non- 
linear bearings, b) backlash, linear bearings, c) backlash, non-linear 
bearings; corresponding linear response. 
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simultaneously in the dynamic model. 


3.4.5. Internal versus External Excitation 

Here we will apply only one type of excitation at a time and compare frequency 
response characteristics of the system due to sinusoidal excitations (F(t)}j and {F(t)) e 
in equation (3.3a). Here, only the bearing non-linearity will be considered assuming 
no tooth separation, 1^=0. In the case of internal excitation, amplitude of alternating 
component is dependent on 12 h . Conversely, the amplitude of alternating external force 
is frequency independent. Therefore, for a constant force ratio F ai /F m , i=h,T, the 
overall amplitude ratio f° r internal excitation varies with excitation frequency 

whereas the overall amplitude ratio F a j/F m for external excitation remains constant 
similar to the single degree of freedom model of Chapter II. 

Figures 3.10a and 3.10b compare y g 2 a versus Cl and p a versus 12 spectra for 
12 t =12 h =12,F m =1.0, F aT =F ah =0.5, F b2 =0, C 22 =0.01, C 23 =0.0125, C 33 =0.05 and 
K 22 = K 23=0.25. As shown Figures 3.10a and 3.10b, both excitations yield the same 
values at 12=1.0 since F m /F aT =F m /Q 2 F ah at 12=1.0. For 12>1 the internal excitation 
gives larger amplitudes than the external excitation, but the converse is true for 12<1. 
Another important difference is that for increasing 12, there are two jumps ( a jump-up 
from the no-impact regime to the double sided impact regime at 12=0.3 and a jump- 
down from the double sided to the single sided impact regimes at 12=0.4) for the 
external excitation; but the double sided impact regime is not seen in the internal 
excitation case and only one jump from the no-impact to the single sided impact is 
found at 12=0.35. 
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frequency 


Figure 3.10. Comparison of frequency responses due to internal static transmission 
error and external torque excitation; F m — 1.0, F^— 0.5, F|j 2-0, bh-0, 

b b2= b c< C33= 0 05 > C23=° ° 125 - C22=°- 01 > k = 1; a) y g2a versus Q h , b) 
p a versus Qjj. 
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3.5. THREE DEGREE OF FREEDOM SYSTEM STUDIES 

3.5.1. Classification of Steady State Solutions 

The steady state response results excited by internal force {F(t))j at each 
frequency have been categorized into the following five groups: (i) harmonic or 
nearly harmonic solution at the same period t p =27t/£2| 1 as that of forcing function 
(period-one), (ii) non-harmonic period-one solution, (iii) subharmonic solution with 
period ntp, n>l, (iv) quasi-periodic solution, and (v) chaotic solution (non-periodic, 
n — > oo). The solution classification criteria are based on time histories, phase plane 
plots, Poincare maps and Fourier spectra [66-71]. Figures 3.11 through 3.15 illustrate 
different types of steady state solutions obtained from the three degree of freedom non- 
linear model of Figure 3.2a. As shown in Figure 3.11, period-one, non-harmonic 
solutions have a non-circular (non-elliptic) phase plane plot and repeat themselves at tp. 
Figure 3.12 shows the time histories, phase plane plots and Poincare maps for a 2t p 
subharmonic solution. In this case, the Poincare map consists of two discrete points. 
Similarly, the 6t p (ultra)subharmonic response has six points on the Poincare map as 
shown in Figure 3.13. Figure 3.14 illustrates a typical quasi-periodic response 
(combination oscillations) which consists of two or more "incommensurate 
frequencies [66,67]. Quasi-periodic solutions result in closed orbits on the Poincare 
map as shown in Figure 3.14c. The chaotic responses which are defined by a non- 
periodic time history and as many points as the number of cycles considered in the 
analysis on the Poincare map are illustrated in Figure 3.15. Figures 3.16a-e show the 
FFT spectra corresponding to the time histories given in Figures 3.1 1 to 3.15. For 
period-one non-harmonic solution of Figure 3.11, corresponding spectrum has peaks at 
mco where to is the fundamental frequency and m is a positive integer, as shown in 
Figure 3.16a. The nt p subharmonic solutions contain peaks at the frequencies mco/n. 
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For instance, spectra of the 2tp and 6tp solutions of Figures 3.12 and 3.13 have peaks 
at frequencies mco/2 and mco/6 respectively in Figures 3.16b and 3.17a. The quasi- 
periodic solution of Figure 3.14 consists of two fundamental frequencies (Oj and <x >2 at 
a ratio C 0 j /to 2 =8 (an irrational number) and there are peaks at the combination 
frequencies mt 0 j+rC 02 , m,r = 0, ±1,±2 — Finally, a characteristic broad band 

spectrum is obtained in Figure 3.18 when the solution as shown in Figure 3.15 is 
chaotic. 

3.5.2. Routes to Chaos 

In Chapter n, we had investigated the effect of mean load F m and damping £ on 
chaotic response of a gear pair. It has been shown that chaos typically exists in 
lightly damped and lightly loaded gear pairs. Now, two different routes to chaos for 
the three degree of freedom system are illustrated. 

A: Perjqd- Doubling Route to Chaos : This consists of a sequence of bifurcations 
of the periodic response to another periodic response with twice the period of the 
original response due to a change in one of the system parameters [66, 68-70]. Figures 
3.19 to 3.23 demonstrates this with a change in excitation frequency ft h given F m =0.1, 
Fah = 0-05, F bi =0 i— 1,2, ^11-^22=0 01, C 13 =C 2 3=0.0125, C 33 =0.05, K 11 =k 2 2=1 25, 
k 13= k 23=°-25, b c =b h 311(1 b bi=°. 1=1 >2- At ft h = 1.5 00, with zero initial conditions, 
we note a t p harmonic solution in Figure 3.19. This tp solution bifurcates to 2t p 
subharmonic solution in Figure 3.20 when ft h is reduced to ft h =l .48. Furthermore at 
ft h = 1.44 a 4t p solution is obtained which is further transformed to a 8tp response at 
ft h =1.402 as shown in Figures 3.21 and 3.22 respectively. At ft h =1.4, however the 
steady state solution becomes chaotic as evident from Figure 3.23. Also note from 
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Figure 3.18. FFT spectra of yg 2 (t) for the chaotic response given in Figure 3.15 
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spectra of Figures 3.19 to 3.23 that the single peak spectrum is transformed to a broad 
band spectrum corresponding to chaos through period-doubling. 

B. Ouasi-periodic Route to Chaos: In this case, instead of going through the 
period doubling, the response undergoes a sequence of Hopf bifurcations with a 
change in a system parameter. First a quasi-periodic response and finally a chaotic 
response are obtained [66,71]; Figures 3.24-3.27 shows such a route to chaos. 3tp 
solution at £2 b = 1 . 1 bifurcates to a quasi-periodic solution when Q b is increased to 1.2 
as shown in Figure 3.24 and 3.25. A further increase in ft b to 1.3 yields a deformed 
closed orbit Poincare map as shown in Figure 3.26. It is then transformed to a chaotic 
strange attractor at £2 b =1.4 in Figure 3.27. Similar observations are evident from the 
Fourier spectra shown in Figures 3.25-3.27; the spectrum consists of combination 
frequencies for the quasi-periodic solution, but it finally changes to a broad band 
spectrum characteristic of chaos. 

3.5.3. Parametric Studies 

The three degree of freedom non-linear model of Figure 3.2a is used to study the 
effect of several system parameters such as the stiffness ratio, k = k bi/ k h i = l>2, 

radial bearing preload to mean load ratio, F b j/F m (i=l,2) and alternating load to mean 
load ratio FjQj/F,,, on the non-linear dynamic behavior. A geared rotor-bearing system 
of one-to-one gear ratio (v g =l) with mgj=l kg, Igj=0.0008 kg/m 2 , dgi=0.08 m, i= 1 ,2 , 
and mesh stiffness k b =2xl0® N/m is selected. Five percent mesh damping and one 
percent bearing damping values, which are somewhat realistic [72,73], are used, i.e. 
^jj=^ 22=0 01. ^i 3=^23 = ® 0125 and £33=0. 05. Both non-linearities are considered 
separately in this parametric study. The case when both bearing non-linearities and 
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Figure 3.19. Frequency doubling route to chaos; F m =0.1, F a i,=0.05, Fj, 2=0, bh=b c , 
b b2 =0 , ^33=0.05, ^3=^23=0.0125, ^11=^22=0 01, k=5; period-one 
response at t2 h =1.5. 
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Figure 3.20. Frequency doubling route to chaos; F m =0.1, F a j,=0.05, F b 2=0, b b =b c , 

bb 2 = 0 , ^33=0.05, ^13=^23=0.0125, Cl l = ^22 = ® -01, k=5; period-two 
response it £2 b =1.48. 
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Figure 3.21. Frequency doubling route to chaos; F m =0.1, F ah =0.05, F b2 =0, b h =b c , 

b b2=°, C33= 0 05 ’ Cl3=C23 =00125 ’ Cl 1=C22= 001 > ^=5; period-four 
response at £2^=1. 44. 
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Figure 3.22. Frequency doubling route to chaos; F m =0.1, F ah =0.05, F b2 -0, b h -b c , 
^2=0^33=0.05, ^13=^23=0.0125, ^11=^22=0.01, k=5; period-eight 
response at f2 b = 1 .402 . 
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Figure 3.23. Frequency doubling route to chaos; F m =0.1, F ah =0.05, F b2 =0, b h =b c< 

bb2—0 ^ 33 — 0.05, ^13=^23=0.0125, Cl 1 = C22 = 0 01, k=5; chaos at 
Q h =1.4. 
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Figure 3.24. Quasi-periodic route to chaos; F m =0.1, F a (,-0.05, F b 2=0, bh=b c =0.1, 

b b2=°, C33= 0 05 > Cl3=C23= 00125 > Cl 1=C22= 001 . ^=1; period three 
response at £2^= 1 • 1 • 
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Figure 3.26. Quasi-periodic route to chaos; F m =0.1, F ab =0.05. F b 2=0, b b =b c =0.1, 


b b2 =0 ^33=0.05, ^ 3 =^ 23=0 0125, ^ 11 =^ 22=0 01, k=l; transition from 
quasi-periodic to chaotic response at Q b =1.3. 









Figure 3.27. Quasi-periodic route to chaos; F m =0.1, F ah =0.05, F b2 =0, b h =b c =0.1, 

bb2 = 0, C 33 = 0 05, ^i 3 = ^ 23 = 0 0125, Cl J = C 22 = ® ^1. k=l; chaos at 
^ h =1.4. 
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backlash exist simultaneously is not considered here as it has already been examined in 
Section 3.4 using the two degree of freedom model. 


A dear Backlash. Linear Bearings : For the non-linear model of Figure 3.2a with 
linear bearings (b bi =0), bt, is used as the characteristic length b c . Figures 3.28a and 
3.29a show y gla (Q h )and p a (fl h ) spe ctra respectively for the soft bearings (k = 1) with 
F m =0.1 and Fah/Fm=°- 5 - Here y g 2a spectrum is not included since y g i a =y g 2a <y g l (t)= " 
y g2 (t)) for a one-to-one gear ratio. As evident from Figures 3.28a and 3.29a, y g i a at 
the first peak is larger than that seen at the second peak. Conversely, p a is much lower 
at the first peak, implying that first mode is dominated by the transverse vibrations 

whereas the second mode is dictated by the torsional vibrations for soft bearings. This 
corresponds to a large dynamic bearing force F dbi »k bi y gi , but a small dynamic 

transmission error (p(t)+e(t)) around the first natural frequency. Jump phenomenon is 
observed only at the second peak which is governed by the gear mesh where the 
excitation e(t) is applied. This indicates that for k = 1 the modes are weakly coupled. 
Accordingly the gear mesh non-linearity, which forces the second peak to be non- 
linear, does not affect the linear characteristics of the first peak. However the modal 
coupling becomes stronger with an increase in k as evident from the jumps seen at both 
peaks in Figures 3.28b and 3.29b for k = 5 and Figures 3.28c and 3.29c for k = 10. 
Torsional vibrations start to dominate the first mode, and the second peak shifts upward 
and it may eventually move out of the operational speed range. Therefore, for a large k 
value associated with stiff shafts and bearings, a single degree of freedom torsional 
model as used in Chapter II should be sufficient. Another effect of a large k is that the 
nt p subharmonic solutions replace some of the harmonic solutions as shown in Figures 

3.28b,c and 3.29b,c. 
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Figure 3.28. y gia versus^, plots for F m =0.1, F ah =0.05, F bi =0, b h =b c ; 

bbi _ ®» a) k=l, b) k=5 , c) k=10; a harmonic, A period-one 
harmonic, o period-n subharmonic, • quasi-periodic or chaos, — 
corresponding linear response. 
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Figure 3.29. 



g a versus plots/or F m =0.1, Fah-0.05, Fy— 0, bj,— b c -0.1, b^i 0, a) 
k=l, b) k=5, c) k=10; <=> harmonic, A period-one non-harmonic, o 

period-n subharmonic, • quasi-periodic or chaos, corresponding 

linear response. 
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Figure 3.30. y™ versus Q h plots for F m =0.1, F bi =0, b h =b c =0.1, b bi =0 £ 33 =0.05, 

Cl3=C23= 0012 5. Cl 1=^22 =001. k =l; a) F ah =0.05, b)' F ah =1.0; 

harmonic, o period-n subharmonic, • quasi-periodic or chaos, 

corresponding linear response. 
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Figure 3.31. 


p a versus Q b plots for F m =0.1, F ah =0.05, F b j=0, b h =b c =0.1, 1^=0, a) 
k=l, b) k=5, c) k=10; ° harmonic, o period-n subharmonic, • quasi- 
periodic or chaos, corresponding linear response. 



Figures 3.30a, 3.31a and 3.30b, 3.31b examine the effect of load ratio for 

A 

F a h/F m =0.5 and 1 respectively with F m =0.1 and k = l. The frequency response is 
well-defined and dominated by the harmonic solutions for F a j 1 /F m =0.5. However 
when Fjj/Fjn is increased to 1, the region beyond flj,=1.0 becomes strongly non-linear 
consisting of subharmonic, quasi-periodic or chaotic responses. This was also seen for 
the gear pair model. 

B, No Gear Backlash. N on - linear Bearings: Now consider the same system with 
bh=0 and non-linear bearings described by the approximate model given in Sections 
3.2.4. Results are presented for both roller and ball bearings separately since the radial 

clearance in a ball bearing is much smaller than that typically found in a roller bearing 
with the same inner and outer race diameters, say bn; * lObw; • For each case, 

J “roller “ball 

bj,i is used as the characteristic length b c . 

Figures 3.32 and 3.33 show the steady state frequency response plots for the 
roller bearings of b b i=0.01 mm (i=l ,2), for three different k = k bi / k h values. As 

shown in Figures 3.32a and 3.33a for soft bearings (k = 1), two modes of interest are 
weakly coupled, and therefore the jump phenomenon is seen only at the first peak 
which is dictated by mostly transverse vibrations. Here, the bearing non-linearities 
affect the transverse displacements y g j and yg 2 which makes the first mode non-linear. 
The second peak is more dependent on the modal interactions. Here, again all of the 
solutions are purely harmonic and only no impact and single sided impact regimes 

A 

exist. However, an increase in k introduces chaotic and subharmonic responses as 
shown in Figures 3.32b,c and 3.33b,c. For instance, all the solutions within the range 
0.6<X2jj<l .7 are non-harmonic for k = 10; the frequency response is no longer well 
defined and modal coupling is sufficiently strong so that jump discontinuities are 
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observed at both peaks. 

A decrease in Fm/Fjj, enhances the degree of the non-linearity as shown in 
Figures 3.34a, b and 3.35a, b. At Fm/F^l in Figures 3.34b and 3.35b, the responses 
within 0.5<£2jj< 1.4 are scattered and non-harmonic. However radial preloads F bi on 
bearings should help the mean load F m in limiting the effect of non-linearities. To 
prove this claim, we apply a high radial preload to both bearings, say F m =Fbi,i = U, 
corresponding to F b j=0.25 for F m = 1 .0 in the dimensionless form, and compare results 
with the previous case of F bi =0 in Figures 3.34a, c and 3.35a, c. Figures 3.34c and 
3.26c show that most of the subharmonic and chaotic solutions of Figures 3.34a and 
3.35a are replaced by harmonic solutions. Hence a well defined frequency response 
curve with clear jump discontinuities is obtained by applying F bi . 

Now, replace roller bearings with ball bearings with b c =bbj=0.001 mm and again 
investigate the effect of k = k w /k h in Figures 3.36 and 3.37. With soft bearings 

(k = 1) the frequency response is linear (no jumps) which indicates that there are no 
impacts within the frequency range considered as shown in Figures 3.36a and 3.37a. 
However at a larger value of k , double sided impact solutions appear since the 
clearance 2b b j is very small. In Figures 3.36b and 3.37b for k = 5, the frequency 

A 

region of 0.65<f2 b <0.8 consists of the double sided solutions. When k is increased 
to 10, the same behavior is also seen at the second peak as shown in Figure 3.36c and 
3.37c. The earlier discussion on the effect of Fab/Fm and F bi /F m for roller bearings is 
also applicable to the ball bearings as well. 

3.6. EXPERIMENTAL VALIDATION 

Munro's experimental dynamic transmission error results [27], which were 
acquired in 1962 using a four-square spur gear test rig, have been used earlier to 
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Figure 3.32. y g j a versus plots for the case of rollerbearings; F m =1.0. F^O.5, 
Fbi=0, bb=0, bbi=b c =0.01mm, a) k=l, b) k=5, c) k=10; <=> harmonic, A 
period-one nonharmonic, o period-n subharmonic, • quasi-periodic or 
chaos, corresponding linear response. 
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Figure 3.33. p a versus Q h plots for the case of roller bearings; F m =1.0, F ah =0.5, 
Fbi=0, bh=0, bbi=b c =0.01mm, a) k=l, b) k=5, c) k=10; a harmonic, A 
period-one nonharmonic, o period-n subharmonic, • quasi-periodic or 
chaos, corresponding linear response. 





Figure 3.34. ygj a versus O h plots for the case of roller bearings; F m = 1.0, bj,=0, 
b b i=b c =0.01mm, k=5; a) F ah =0.5, F bi =0, b) F ah =0.5, F bi =0.25, c) 
F ab =10, F bi =0; a harmonic, A period-one nonharmonic, o period-n 

subharmonic, • quasi-periodic or chaos, corresponding linear 

response. 
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Figure 3.35. p a versus ft b plots for the case of roller bearings; F m =1.0, b h =0, 
b bi =b c =0.01mm, k=5; a) F ab =0.5, F b j=0, b) F ab =0.5, F b p0.25, c) 
F ah =1.0, F bi =0; a harmonic, A period-one nonharmonic, o period-n 

subharmonic, • quasi-periodic or chaos, corresponding linear 

response. 



Figure 3.36. ygj a versus £2}, plots for the case of deep-groove ball bearings; 

F m =10, F ah =5, F bi =0, b b =0, b bi =b c =0. 001mm, a)k=l, b)k=5, c) 
k=10; a harmonic, A period-one nonharmonic, o period -n subharmonic, 
• quasi-periodic or chaos, corresponding linear response. 
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Figure 3.37. p a versus fl h P lots for the case of deep-groove ball bearings; F ra -10, 
F ah =5, F bi =0, b h =0, b bi =b c =0. 001mm, a) k=l, b) k=5, c) k=10; 
harmonic, A period-one nonharmonic, o period-n subharmonic, • quasi- 
periodic or chaos, corresponding linear response. 
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validate the spur gear pair model of Chapter II assuming that the modes are weakly 
coupled with each other. Now, we use the three degree of freedom non-linear model of 
Figure 3.2a to analyze his test set-up [27]. In our analysis, bearings are assumed to be 
linear since they were highly preloaded; only gear backlash non-linearity is considered 
with 2b|j=0.12 mm [74]. Table 3.2 lists the system parameters used in equation (3.3a) 
under four different mean loads. The system damping ratios are assumed to be uniform 
at each load. Average mesh stiffness kj, and alternating load values associated with 
each mean load F m arc also tabulated in Table 3.2. 

Figures 3.38 through 3.42 compare dynamic transmission error predictions with 
the measurements [27]. At the design load (DL), which corresponds to the minimum 
static transmission error, good agreement is seen except for the jump discontinuity 
found experimentally at the second peak. In our model, we have increased slightly 
beyond the value given in Reference [27], so that the predicted amplitude of the first 
linear peak matches very well with the experiment. Such slight changes in critical 
system parameters such as F^ and although they are varied within the experimental 
uncertainties, may alter the frequency response drastically as illustrated in Figures 
3.39a and 3.39b. In Figure 3.39a we note that when we reduce the force ratio Fj^/F^ 
from 30 (the value given in Reference [27]) to 10, a jump discontinuity at the second 
peak is seen which is compatible with the experiment. Similarly, a small change in the 
damping values, which are not reported in Reference [27], also affects the frequency 
response as shown in Figure 3.39b. Figures 3.40 to 3.42 compare the dynamic 
transmission error spectra at 3/4, 1/2 and 1/4 of the design load. From these spectra we 
conclude that our proposed theory agrees well with the experimental results of Munro 
[27] both qualitatively and quantitatively, although amplitudes in the off-resonance 
regions are slightly off. Such amplitudes are close to the noise floor in the experiment. 
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Table 3.2. Parameters of the Munro's experimental set-up extracted from Reference 


[27]. 




Mean load, F m 


design load(DL) 
0.183 

3/4 of DL 
0.146 

1/2 of DL 
0.105 

1/4 of DL 
0.0579 

Fah 

0.0058t 

0.0178 

0.0296 

0.0393 

k b (N/m) 

1.16xl0 9 

1.16x109 

1.16x109 

1.16x109 

kjj (N/m) 

3.44xl0 8 

3.22xl0 8 

3.01xl0 8 

2.72xl0 8 

K ll» k 22 

0.950 

0.966 

0.983 

1.007 

k 13* k 23 

0.242 

0.242 

0.242 

0.242 

Cll> *»22^ 

0.01 

0.01 

0.01 

0.01 

Cl 3’ ^23 t+ 

0.00375 

0.00375 

0.00375 

0.00375 


0.015 

0.015 

0.015 

0.015 


t modified so that the linear peak matches with the measured value, 
tt estimated using the experimental data of Reference [27]. 
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Figure 3.38. Comparison of theory with Munro's [27] experiment at the design load. 
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Figure 3.39. Effect of a) alternating load and b) damping on the dynamic transmission 
error spectrum at the design load. 
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Figure 3.42. Comparison of theory with Munro’s [27] experiment at 1/4 of the design 
load. 
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Further, discrepancies between the theory and the experiment may be due to the 
assumptions made in this study. It should however be pointed out that the experiment, 
even though it was very precise, was conducted almost three decades back before the 
advent of modem non-linear dynamics and chaos science [66-71]. Interestingly, 
Munro [27] had also reported subharmonic responses at 3/4, 1/2 and 1/4 of the design 
load, and "unrepeatable responses" (probably chaos) at 1/4 of the design load. Our 
results agree with such experimental observations as well. Both experiment and theory 
show that the chaotic solutions dictate the frequency response beyond the second 
resonance for 1/4 of the design load in Figure 3.42. 

3.7. CONCLUSION 

This analytical study on the non-linear dynamics of a geared rotor-bearing system 
with gear backlash and bearing clearances, as excited by the internal static transmission 
error and/or external torque pulsations, has made a number of contributions to the state 
of the art. First, an approximate non-linear model of the rolling element bearings with 
clearances has been developed and validated. Second, our mathematical model has 
been validated as it compares well with a previous experimental study, and several key 
parameters such as the mean load, mean to alternating force ratio, radial bearing 
preload, bearing stiffness and damping have been identified. Third, conditions for the 
chaotic and subharmonic resonances, and the routes to chaos have been identified. 
Forth, non-linear modal interactions have been examined. Finally, on a more 
fundamental note, our study enriches the current literature on coupled vibro-impact 
pairs [52,54,55]. 
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CHAPTER IV 


INTERACTIONS BETWEEN TIME-VARYING MESH STIFFNESS 
AND CLEARANCE NON-LINEARITIES 

4.1. INTRODUCTION 

Dynamic models of geared systems can be classified into four main groups. The 
first group includes linear time-invariant (LTI) models, as evident from an extensive 
review of the literature given in Reference [14]. The second group considers linear 
time-varying (LTV) mesh stiffness k h (t) in the analysis [27,29,62]. Periodic variation 
in kj,(t) is due to the changing of the number of conjugate teeth pairs in contact during 
the convolute action. Accordingly, the system is excited parametrically as well as by 
the static transmission error e(t) introduced by kinematic errors and tooth deflections. 

In this case, the equation of motion of the gear pair essentially reduces to the Mathieu's 
or Hill’s equation with a periodic external forcing function. The third group includes 
gear backlash in the models, but with time-invariant average mesh stiffness k^ * kj,(t) 

[25,30,31,33,34,41,42]. It should be noted that backlash is bound to exist either by 
design or due to manufacturing errors and/or wear in any gear pair. Finally, the last 
group considers both gear backlash and mesh stiffness variation simultaneously 
[24,35-37,53]. However, none of these studies have addressed explicitly the effect of 
k h (t), including its interaction with the backlash non-linearity, on the steady state 
frequency response. Ozguven and Houser [24] have attempted to analyze this problem 
by replacing k h (t) with a constant mesh stiffness and by defining the "loaded static 
transmission error" excitation at the mesh point. But, Ozguven [75] has stated recently 
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that this approach may not work depending on the system parameters, and 
recommended a detailed investigation of this issue. 

Numerous publications have analyzed LTV systems [76-78] and time-varying 
systems with quadratic and/or cubic non-linearities [79-82]. However, such studies are 
not directly applicable to the geared rotor systems which are the main focus of this 
chapter. This problem requires the solution to a set of time- varying differential equations 
with clearance-type non-linearities as excited by a periodic force generated at the gear 
mesh. 

4.2. PROBLEM FORMULATION 

A generic geared rotor-bearing system, which consists of a spur gear pair 
mounted on flexible shafts which are supported by rolling element bearings as shown 
in Figure 4.1a, is considered here. The gear box is assumed to be rigid. The effect of 
the prime mover or the load inertia is not considered assuming that such inertial 
elements are connected to the gear box through soft torsional couplings. Further, we 
assume that the system is symmetric about the plane of the gears and the axial motion 
(parallel to the shafts) is negligible. The governing equations of motion can be given in 
the matrix form as 

[M]{q"(t)} + [C]{q'(t)} + [K(t)]{f(q(t))} = {F(t)} (4.1) 

where [m] is the time -invariant mass matrix and {q(t)} is the displacement vector. 
Here, damping matrix [c] is assumed to be time-invariant, as the effect of the tooth 

separation and time-varying mesh properties on mesh damping are considered 
negligible; validity of this assumption will be examined later. The stiffness matrix 
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(a) 



(b) 


a) Generic geared rotor-bearing system, b) clearance non-linearity in 
gears and bearings. 
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[K(t)] is considered to be a periodically time-varying matrix given by 
[K(t)] = [K(t + 2rc/Qh)] where ft h is the fundamental gear mesh frequency. The 

non-linear displacement vector {f(q(t))j includes the radial clearances in bearings and 
the gear backlash as shown in Figure 4.1b, and the forcing vector {F(t)} consists of 

both external torque and internal static transmission error excitations. 

This chapter extends our previous non-linear single degree of freedom spur gear 
pair model of Chapter II and multi-degree of freedom geared rotor-bearing system 
model of Chapter ID by including time-varying mesh stiffness k h (t), and investigates 
its effect on the frequency response of lightly and heavily loaded geared systems. 
Interactions between mesh stiffness variation and system non-linearities associated with 
gear backlash and radial clearances in rolling element bearings are also considered. 
Resonances of the corresponding LTV system associated with the parametric and 
external excitations are identified using the method of multiple scales. Our formulation 
will be validated by comparing predictions with available experimental results [27,28]. 

4.3. MATHEMATICAL MODEL 

A reduced order form of the multi-degree of freedom system of equation (4.1), 
the three degree of freedom non-linear model of the geared rotor system, which has 
been used in Chapter IH, is considered here. This includes gear inertias I g i and I g2 , 
gear masses m gl and m^, and base circle diameters d gl and d g2 , as shown in Figure 
4.2a. The gear mesh is described by a non-linear displacement function f b with time 
varying stiffness k h (t) and linear viscous damping c h . Bearings and the shafts that 
support the gears are modeled by equivalent elements with viscous damping 
coefficients c bl and c b2 and non-linear springs defined by force-displacement functions 
f bl md f b2 which are approximated by piece-wise linear, dead zone type non-linearities 
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as suggested in Chapter ID. Both low frequency external excitation due to torque 
fluctuations and high frequency internal excitation due to the static transmission error 
e(t) are considered in the formulation. Input torque fluctuation is included, but the 
output torque is assumed to be constant, i.e. T gl (t) = T glm + T gla (t) and 

T g2 (t) = T g2m . External radial preloads and fy 2 are also applied to both rolling 
element bearings. 


4.3.1. Equations of Motion 

Equations of coupled transverse-torsional motion of the geared rotor-bearing 
system shown in Figure 4.2a with the displacement vector 
{q(0} = {y g i(t), y g2 (t), p(t)} are given in matrix form as follows: 


" m gl 

0 

0 

ygi <*>' 


’ c bi 

0 

c h 

y g i(t) 

0 

m g2 

0 

yg2(0 

► + 

0 

c b2 

-c h 

yg2(0 • 

“ m cl 

m cl 

m cl. 

. p"(t). 


0 

0 

c h. 

v'd) . 


k bl 0 k h (t)' 

f bi(y g i(0) 


-*bl 

0 k b2 -k b (t) 

f b2(y g 2(0) 

► — * 

hi 

0 0 k h (t)_ 

. f h (p(t)) . 


F m -m cI e"(t) + F aT (t) 


(4.2a) 


p( t ) = 0 g i (t ) - -jp- 0 g2 (t ) + y g i ( t ) - y g2 (t) - e(t ) ; 


(4.2b) 


mcl * / .2 


1 


A t A 2 A 

d 8l , d g2 


4 ! g i 4 I g2 J 


■ p _ 2T glm _ 2T g 2m = _ m cl T gla( t ) 

’ r m ~ , l aTV t ) = ; ( 4 . 

d gl d g2 2I gl 


2c-e) 
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Figure 4.2. Non-linear models of the geared system of Figure 1; a) three degree of 
freedom system, b) gear pair model with a single degree of freedom. 
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(4.20 


^ oo 

k h (t) = k h (t + =-) = k hm + Ikhaj cosO-Qfaf-H^); 
w h r=l 

Jygi-t’bi; y p > b bl f(p-b„); p>b„ 

fbi(ygi>-|0; - b b , < ygi < b 6i ; f(,(p)— • 0 ; -b h <p<bb (4.2g.h) 

ygi + ^bii ygi < — ^bi (P + *> h ); p<-b h 

Here, ( )' donates derivative with respect to time t , y gi and 0^ are the transverse and 
torsional displacements of the i-th gear (i=l,2), m*., is the equivalent gear pair mass, 
F m is the average force transmitted through the gear mesh, and F aX (t) is the fluctuating 
force related to the external input torque excitation. Equations (4.2a-h) have been used 
previously in Chapter III except for kj, (t) term, which is expanded in equations (4.2a- 
h) in the Fourier series form. Here p(t) is the difference between the dynamic 
transmission error and the static transmission error e(t). A dimensionless form of 
equation (4.2a) is obtained by letting yg|(t) = yg|(t) / b c , p(t) = p(t)/b c , 

®n — V^hm ~i m cl > w bi — sj^bi / m gi (i = l>2) and t = to n t where b c is the characteristic 
length. Here, we consider periodic excitation for both e(t) and F aT (t) as 

po oo 

®(*) = X.e r c °s(rQ h t + 4> er ), F aT (t)= Z F a Tr cos(rfl T t + <t>Tr ) where and fXr are 
r=l r=l 

the fundamental excitation frequencies of internal displacement and external torque 
fluctuations, respectively. Further, define dimensionless excitation frequencies 

to yield the following governing equations of 

motion in the dimensionless form: 
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-Fbl 


0 


+ 

0 

k 22 

“^23(0 

fb2(y g 2> 

► zzz * 

Fb2 

> + < 

0 

>+' 


0 

0 

k 33(0. 

. f b(p) . 




FahW 



Fah(t)= I FahrCrfih^cosCrQht + ilJer); 
r=l 



00 r F a Tr 

F aT (t) = I F aTr cos(rf2 T t + <t> Tr ); F aTr = “T ; 

r=i m cl b c®n 


F m = 


L m 


; F b. = -r^ !i TT' i = 1 ’ 2 ' 


/ 9 1 ^ 

m gi b c w n 


Cti = Cbi ; Ci3= C - — , i=l,2; ^ 33 = Cfa 

2mgi0) n 2m^03 n 2m cl (D n 


©hi „ _ ^hO) • 1 *j, 
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K 33( t ) = 7^ = l+ XerCO^r^ + ^hr); 
^hm r=l 


0 

0 •; (4.3a) 
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(4.3b, c) 


(4.3d,e) 
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(4.3m) 
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where F bi (i=l,2) and F m are the dimensionless components of the mean force vector 
{Fj m , and F aT (t) and F^t) pertain to the alternating external excitation {F(t)} e and 
internal excitation {F(t)}j force vectors, respectively. 

4.4. GEAR PAIR STUDIES 

First, we consider the spur gear pair model of Chapter II, as shown in Figure 
4.2b, and investigate the effect of k h (t) on the steady state frequency response. The 

equation of motion of the gear pair is obtained from equation (4.3a) by substituting 
ygi(t)=0 for fixed gear centers. Neglecting the input torque variations for the sake of 
convenience, F aX (t)=0, we obtain the following equation of motion of a loaded gear 

pair with time-varying mesh stiffness and backlash, as excited by the static 
transmission error. 


p(t) + 2^33p(t) + K 3 3(t)f h (p(t)) = F m +F ah (t) (4 4) 

where £33, *33(1), fh(p), F m and F^ are still given by equation (4.3). Since none of 
the existing analytical solution methods are found to be suitable for this problem, the 
digital simulation technique is used to solve equation (4.4). A 5th-6th order Runge- 
Kutta numerical integration algorithm with variable time step [50] is used here. This 
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technique has already been employed successfully in Chapter II for the gear pair 

problem with time-invariant mesh stiffness. 

First consider only sinusoidally varying mesh stiffness and the static transmission 
error excitation, i.e. K 330 ) = 1 + £j cos(ftj,t + 4>bi ) where 6=£j and 
F ab (t) = Fjhjflh 2 cos(Qfjt + <|>ei ) ■ The tooth deflection or F^t) is maximum under the 
applied mean load F m when K 33 (t) is minimum (i.e. minimum number of gear pair in 
contact). Similarly, minimum F^t) corresponds to maximum K 33 (t). Therefore, there 
is an out of phase relationship between F^ft) and K 3 3 (t), i.e. <|>hl =< t>el +7t - Here, we 
set <t>bi=rc and <{> el =0 for the sake of convenience. Hence, equation (4.4) is modified 

to: 


p(t) + 2^33p(t) + [1 — £cos((2ht)]fjj(p(t)) — F m +F a j 1 (2| 1 cos(ftt,t). (4.5) 


4.4.1. Linear Time-varying (LTV) System 

For zero gear backlash b h , the gear mesh displacement function is f h (p)=p. Hence 

equation (4.5) reduces to an LTV equation as follows: 


p(t) + 2^33f>(t) + [1 — £cos((2jjt)]p(t) — F m + Fjijtijj cos((2jjt ) (4.6) 

Approximate analytical solutions of similar LTV differential equations with parametric 
and external excitations are already given in the literature [76,77,79-82]. Therefore, we 
will not attempt to solve equation (4.6) completely. Instead we identify the 
corresponding resonances using the method of multiple scales [76]. A first order 
uniform solution is given by an expansion having the following form where the scalar 
parameter £ is assumed to be very small. 
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p(t;£) = p 0 (T 0 ,Ti) + £p 1 (T 0 ,T]) + 0(e 2 ) 


(4.7) 


where T n — £ t is the n-th time scale. Expansions for the derivatives with respect to t 
is obtained in terms of the partial derivatives with respect to time scales D n = 9/ 9T n 

d / dt = D 0 + £D, + 0(£ 2 ) ; d 2 / dt 2 = D 2 + 2£D 0 Dj + 0(£ 2 ) (4.8a,b) 


Substituting equations (7) and (8) into equation (6) and equating like powers of £, with 
the external force being applied at 0(£°), one obtains 

^oPo Po = F m ■*" Fahl^h cos (^hTo) (4.9a) 

D oPl + Pi = -2D 0 D!p 0 - pD o p 0 + Po cos(fl h T 0 ) (4.9b) 

where ep=2£ 33 . The solution of equation (9a) is given in the complex domain as 

Po = A ( T l )e ,T ° + F m + Ae ,£ * hT ° + cc ; A = / 2(1 - ill ) (4.1 Oa,b) 

where cc represents the complex conjugate terms and i = V^l . At fl h * 1 , we note that 
(1 - f2 h ) will make the amplitude of the response boundless; this is the primary 

resonance since it appears in the first order. Substituting equation (4.10) into equation 
(4.9b) we get 
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DoPi + P! = -2iD 1 Ae lT ° - i|oAe lT ° - i|iAf2 h e lQhT ° +^Ae l(l+£lh)T ° 

+i-Ae i(,_£lh)T ° +-Ae 2,nhT ° +-F m e‘ n,,To +A + cc (4.11) 
2 2 2 

Besides the primary resonance at £2 h = 1, the particular solution of equation (4.1 1) has 
other secular terms when 2£2 h * 1 and £2 h * 2. At 2£2 h - 1, summation of the 
external excitation frequency and the parametric excitation frequency is close to the 
dimensionless natural frequency which is unity. One could also observe resonances at 
£2 h = n, n>2 when higher scales are considered in equation (4.7). In summary, 
equation (4.6) has resonances at £2^ =0.5 and n, where n=l,2,3 — 

Now, we solve equation (4.6) using digital simulation. Figure 4.3 shows the 
steady state frequency response curves p a (£2jj) and p m (0| 1 ) for a lightly loaded system 
with F m =0.1, Fab^O.05, £ 33 =0.05, 1^=0 and four different e values. Note that e=0 
represents the LTI system. For e>0, we do not observe any multi-valued regions and 
jump phenomenon, similar to those seen for the non-linear systems. Mesh stiffness 
variation £ has a negligible effect on the natural frequency which corresponds to the 
largest peak in Figure 4.3a; but an increase £ in amplifies p a in both resonance and off- 
resonance regions. Here p m is no longer uncoupled from p a , and it varies substantially 
in the vicinity of primary resonance at £2^= 1.0, and the parametric resonance at 
£2 h =0.5, especially for a large E confirming the analysis given earlier. Note that our 
study is limited to only these two resonances as an investigation of subharmonic 
resonances at £2 h = n , n>l is beyond the scope of this study. 
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(a) 



frequency 


Figure 4.3. Frequency response spectra of a lightly loaded LTV gear pair with 
sinusoidal e(t) and kjj(t), F m =0.1, F^jsO.OS, £ 33 = 0 . 05 , bj,=0 and four 
different £ values; a)p a versus b) p m versus flj,. 
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4.4.2. Non-linear Time-varying System 

Next we consider the time-varying non-linear system with backlash, given by 
equation (4.5). Figure 4.4 shows p a (Qjj) an< ^ Pm^h) spectra with F m =0.1, F a i=0.05, 
^33=0.05, and four E values. In this case, we notice a jump discontinuity at the 

resonant frequency and a dual-valued region bounded by jump-up and jump-down 
transition frequencies. For E=0, p m * Pm(^b) evident in the no-impact regime, 

similar to the results given in Chapter n. But when E>0, the transition frequencies 
which define the jump phenomenon become smaller, and p a at the jump-down 
frequency grows with increasing E. This indicates that k|,(t) enhances the degree of 

non-linearity associated with the gear backlash. Similar to the results of LTV system in 
Figure 4.3, we again observe in Figure 4.4a a parametric resonance at fl h =0.5, which 
is strong for E=0.2 and 0.4 curves. Figure 4.5 compares the time histories p(t) for 
each £ values at the parametric resonance peak, Ci|,=0.5. For the LTI system (£=0), a 
harmonic no-impact type steady state solution exists as illustrated in Figure 4.5a. With 
increasing £, this solution is transformed into a non-harmonic periodic solution with a 
larger peak to peak value as a result of the parametric mesh stiffness excitation, as 
shown in Figures 4.5b-d. 

Now, consider a heavily loaded gear pair with F m =0.1, F al =0.01 and £ 33 =0.05. 
Figures 4.6a and 4.6b show p a (O h ) 31x6 Pm^h) spectra, respectively. In this case, 
tooth separation does not occur for E=0 which results in a linear frequency response 
curve with a constant p m . A small jump is seen at £=0.1, and this jump becomes larger 
for £=0.2 and 0.4; overall alternating amplitudes p a at resonance and within the off- 
resonance regions increase considerably. This suggests that the mesh stiffness 
variation is especially important for a heavily loaded gear pair with backlash. 
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frequency 


(b) 



Figure 4.4. Frequency response spectra of a lightly loaded non-linear gear pair of 
Figure 2bwith sinusoidal e(t) and khO), F m =0.1, Fj^^O.OS, £ 33 = 0 . 05 , 
bt,=b c and four different e values; a)p a versus Of,, b) p m versus 
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(a) 



frequency 


Figure 4 . 6 . Frequency response spectra of a heavily loaded non-linear gear pair with 
sinusoidal e(t) and k h (t), F m = 0 . 1 , FahjsO.Ol, £33=0.05, b h =b c and four 
different e values; a)p a versus b) p m versus 
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4.4.3. Non-linear Time-varying Mesh Damping 

In order to investigate the effect of non-linear, time-varying gear mesh damping 
on the steady state frequency response, we consider a sinusoidally varying gear mesh 
damping £33(0 which is assumed to be proportional to the mesh stiffness function 
jcjjft), and a non-linear velocity function gj, associated with tooth separation. Hence 
the governing equation of motion is 

p(t) + 2£ 33 (t)g h (p(t)) + [1 - ecos(fl h t)]f h (p(t)) = F m + Fahflh cos(Q h t); (4. 12a) 


£33(0 = C33 k 33<0 = C33C 1 + ecos(f2 h t + <t> hl )) ; 


(4.12b) 


gh(p( t » = 


0 ; 

b c F b c 
p(t); else 


(4.12c) 


Figure 4.7 compares frequency response spectra for the cases of a) non-linear time- 
varying gear mesh damping defined by equation (4.12), and b) linear time-invariant 
mesh damping, £33 *£33(1) with gh(P( t » = P(0. 355 8 iven b ? eq uati on (4.5). As 
shown in Figure 4.7, both mesh damping models yield virtually the same spectra. 
Therefore a linear time-invariant gear mesh damping model can be used without losing 
any accuracy. 

4.4.4. Periodic e(t) and k h (t) 

Up to now, we have only considered the sinusodally varying kj^t) and e(t)- In 
real geared systems, however, both k|j(t) and e(t) are periodic which can be 
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(a) 



(b) 



Figure 4.7. Frequency response spectra of gear pair corresponding to non-linear time 
varying and linear viscous damping models; a)p a versus Q h , b) p m versus 
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expressed in the Fourier series form. Figure 4.8 illustrates kj,(t) and e(t) for a low 
contact ratio spur gear pair; these typical time histories are predicted using an existing 
spur gear tooth model [21]. The first six Fourier coefficients of k h (t) and e(t) are 
given in Table 4.1. Note that n-th harmonic of e(t) causes a resonance around 

£2 h =l/n. Hence higher harmonics are important only at low frequencies. Accordingly, 

only the first three Fourier coefficients of e(t) are deemed sufficient here, given the 
nature of the practical gear noise and vibration problems. Similarly, is also 

truncated for the sake of convenience, retaining the first three terms. Higher terms 
(n>3) could be easily included in the analysis if necessary. Therefore the governing 
equation with F a -p(t)=0 is 


3 3 2 

p(t) + 2 £ 33 (t)p(t) + [l- Ie r cos(rQ h t)]f h (p(t)) = F m + X^ahr^h) cos(rf2 h t) 

r=l r=l 


(4.13) 


Figures 4.9a and 4.9b show frequency response spectra for F m =0.1, F a hi~0.05, 
F ah2= 0 02 ’ F ah3= 001 - £33=0.05, e l/ e 2= 2 » 311(1 four different £i=e values. 

As shown in Figure 4.9, two other jump discontinuities, which are clearly associated 
with F ah ? and F a h 3 , are found. These jumps are similar to the fundamental resonance 
peak except they have lower amplitudes. In order for these two jumps to exist, F a h2 
and F ah3 excitations have to be large. Once again, it is evident that the periodic mesh 
stiffness enhances the extent of backlash non-linearity. 
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Figure 4.8 e(t) and k h (t) plots for a low contact ratio spur gear pair [21]. Here, 
number of teeth is 28, diametral pitch is 8, the pressure angle is 20 

degrees, and a 6x1 0‘ 4 in. tip modification starting at 26 degrees roll 
angle is applied. 




Table 4.1 Fourier coefficients of e(t) and k h (t) given in Figure 4.8 


Harmonic n 

e(t) 

(H-m) 

Amplitude 

k h (t) 

(MN/m) 

0 

20.00 

107.8 

1 

3.45 

10.3 

2 

0.26 

5.7 

3 

0.45 

1.3 

4 

0.38 

4.0 

5 

0.27 

1.5 

6 

0.26 

1.4 


i 
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Figure 4.9. Frequency response spectra of a non-linear gear pair with periodic e(t) and 
k h(0, F m =0.1, F ahl =0.05, F ah2 =0.02, F^asO.Ol, £33=0.05, b h =b c , 
e l/ £ 2=2, e!/e 3 =4 and four different e^e values; a)p a versus fl h , b) p m 
versus £2u. 




4.5. GEARED ROTOR-BEARING SYSTEM STUDIES 

Next, we solve the three degree of freedom non-linear model of Figure 4.2a and 
equation (4.3a) using the digital simulation technique for excitation F ah (t) with 

F aT (t)=0. 


4.5.1. Sinusoidal e(t) and k^ft) 

Similar to Section 4.4, assume sinusoidal forms: K 33 (t) = l + ecos(Qht + <t , hl) 
and F ah (t) = F ahl O h 2 cos(O h t + <t> el ). Equation (4.3a) reduces to an LTV matrix 

equation when gear backlash and radial bearing clearances are set to zero, i.e. 
f h (p(t))=p(t) and f b i(ygi(0)=y g i(0. i=l,2. Figure 4.10 shows p a (O h ) and y g i a (O h ) 
spectra of this LTV system under heavily loaded conditions with F m =0.1, F^j— 0.0 1 , 
C 33 =0.05, ; i3 =0.0125 and ^*0.01, i=l, 2. Note that y gla =y g2a since the gear ratio 
is one. As shown in Figure 4.10, frequency response has peaks at flb = ®i = 0-4 and 
ft h =<»ii=1.25 where (Oj and to n are the natural frequencies of the LTI system, 
corresponding to the first two coupled transverse-torsional modes. In Figure 10a, we 
observe a parametric resonance at 2Qjj=(Oj similar to the gear pair model of Section 4.4 
whereas the parametric resonance at 2f2j, = ®II is more obvious in p a (I2jj) spectrum of 
Figure 4.10b for a larger £ value. Again like the LTV gear pair results of Figure 4.3, 
y g i a and p a are amplified with increasing £ in both resonance and off-resonance 

regions. 

Now introduce the gear backlash b(j = b c to this heavily loaded system, 
corresponding results are given in Figure 4.1 1. The most significant effect of k^O) is 
that it interacts with the gear backlash non-linearity to develop a jump discontinuity at 
the second primary resonance peak. This is clearly evident from the fact that either 
backlash or k h (t) alone can not cause a jump for such a heavily loaded system. 
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(a) 




Figure 4.10. Frequency response spectra of a LTV geared rotor-bearing system of 
Figure 2a with sinusoidal e(t) and k h (t), F m =0.1, F^sO-Gl, £ 33 =0.05, 
; 33 =0.05, Ci3=0 0125 and {^=0.01, k u = 0.5, ^=1^=0, i=l, 2, and four 
El=t values; a) ygj a versus £l h , b) p a versus O b . 
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Besides an increase in amplitudes in Figures 4.11a and 4.11b, chaotic and/or quasi- 
periodic solutions, which do not exist in heavily loaded time-invariant systems as 
shown in Chapter ffl, are also predicted in the following regimes: 1.35<fl h <1.65 for 
e=0.2 and 1.25<Q h <1.7 for e=0.4. Similar to the LTV system of Figure 4.10, 
parametric resonances at 2£2jj=G)j and 2Q b =fi>jj are again observed here. 

Next, radial clearances in bearings b bi =b c are considered for a heavily loaded 
system with zero gear backlash, F m =1.0, Fahl 3 ®-!* £ 33 = 0 . 05 , £j3=0.0125 and 
£ii=0.01, i=l, 2. In this case, mesh stiffness variation and non-linearity exist in 
different components. As shown in Figure 4.12, £ again increases the amplitudes, but 
the spectrum shape is essentially the same. This is further evident from Figure 4.13, 
given for a lightly loaded drive with F m =1.0 and Fj^^O.5. This suggests a weak 
interaction between the gear mesh stiffness variation and bearing non-linearities. 

4.5.2. Periodic e(t) and k b (t) 

Similar to the gear pair analysis of Section 4.4.4, we consider periodic K 33 (t) and 
F ab (t) with three Fourier coefficients. Figure 4.14 shows spectra for F m =0.1, 
Fah! = 0.01, F ah2 =0.004, F ah3 =0.002, £ 33 =0.05, £ i3 =0.0125 and £ u =0.01, i=l, 2; 
ei/e 2 =2, £^£2=4, and £^£=0, 0.1 and 0.2. Here only the gear backlash non-linearity 
(b h =b c and 1^=0) is considered. Peaks at fl h =( 0 I /n and li h =con/n, n=l,2,3 are 
predicted corresponding to excitation F^. Periodic mesh stiffness enhances the 
alternating amplitudes over the values given by the time-invariant mesh stiffness case 
and introduces a jump discontinuity at the second primary resonance similar to the 
sinusoidal mesh stiffness case. 


189 



(a) 




Figure 4.11. Frequency response spectra of a non-linear geared rotor-bearing system 
with sinusoidal e(t) and k h (t), F m =0.1, F ahl =0.01, £33=0.05, 
£33=0.05, £i 3 =0.0125 and £ u =0.01, k u =0.5, b,,=b c , b bi =0, i=l, 2, and 
four E|=e values; a) y^ versus ft h , b) p a versus fl h . 
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“ E=0. 

“• e=0.2 



(b) 



Figure 4.12. Frequency response spectra of a non-linear geared rotor-bearing system 
with sinusoidal e(t) and ^(t), F m =l .0, Fjjj^O.l, ^33=0.05, C33— 0-05, 
; i3 =0.0125 and ^=0.01, ^=0.5, bfe^bc, ^=0, i=l, 2, and four e 1= e 
values; a) y^ versus Ci h , b) p a versus ft h . 

191 




(a) 




Figure 4.13. Frequency response spectra of a non-linear geared rotor-bearing system 
with sinusoidal e(t) and k h (t), F m =1.0, F^jsO.S, £ 33 = 0 . 05 , £ 3 3=0.05, 
^3=0.0125 and £^=0.01, KjpO.5, bjjpb,., 1^=0, i=l, 2, and four £j=E 
values; a) ygj a versus b) p a versus Qj,. 
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amplitude amplitude 


(a) 



frequency 

Figure 4.14. Frequency response spectra of a non-linear geared rotor-bearing system 
with periodic e(t) and kh(*). F a hi=0.01, Fah2 = ®-®®^* 

F ah 3=0.002, ^33=0.05, ^33=0.05, ^3=0.0125 and £,1=0.01, ^=0.5, 
b h =b c , b bi =0, i=l, 2, ei/e 2 =2, 6i/e 3 =4 and three e t =e values; a)y gia 
versus ft b , b) p a versus fi b . 
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4.6. EXPERIMENTAL VALIDATION 
4.6.1. Spur Gear Pair Dynamics 

First we compare our theory with experimental results of Kubo [28], as extracted 
from Reference [24]. Kubo used a heavily damped (£ 33 = 0 . 1 ) four-square spur gear 
test rig, and measured dynamic factor as the ratio of the dynamic to static tooth root 
stresses. The experimental set-up was designed to support a gear pair with very stiff 
shafts and bearings. Therefore our gear pair model can represent the test rig 
adequately, similar to Chapter II. The static transmission error e(t) and time-varying 
mesh stiffness k b (t) of the tested gear pair has been predicted using an existing elastic 
spur gear model [21] and then equation (4.4) is solved to predict the dynamic response. 
Here, we define the dynamic factor as the dynamic to static mesh force ratio which is 
equivalent to the dynamic factor calculation based on the stress analysis under the 
assumption that the change in the moment arm due to changes in the contact point is 
negligible. Figure 4.15 compares the envelope of measurements obtained by Kubo for 
several tooth pairs with our predictions. When the mesh stiffness k b (t) is assumed to 
be time-invariant, the predicted jump discontinuity is not as large as the jump seen in 
experimental data and the predicted transition frequency is higher than the measured 
value, as reported in Chapter n. However predictions improve significantly when 
sinusoidally varying k b (t) and e(t) at are used, as given by equation (4.5). A 

sharp jump discontinuity is found which matches well with experiment, and the 
predicted dynamic factor is very close to measured envelope. Finally, prediction agrees 
with experiment even better when periodically varying e(t) and kh(t) with the first three 
Fourier coefficients considered. This figure clearly shows that the time-invariant linear 
or non-linear model can not predict the true dynamic behavior, as time-varying mesh 
stiffness must be included in the non-linear mathematical formulation. 
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Figure 4.15 Comparison of theory with Kubo's [24,28] experimental results. 
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4.6.2. Geared Rotor-Bearing System 

As the second example case, experimental results of Munro [27] are compared 
with our geared rotor-bearing model of Figure 4.2a and equation (4.3). These 
experimental results were used earlier to validate time-invariant multi-degree of freedom 
non-linear system of Chapter m in which we had to reduce the damping ratio and 
increase the excitation e(t) in order to correlate theory with experiment. Now we go 

back to the original system parameters given in References [27] and Chapter III, and 
predict the dynamic transmission error spectra using both time-varying and time- 
invariant gear mesh stiffness formulations for a three degree of freedom non-linear 
model with gear backlash and linear bearings. Figure 4.16 compares results at the 
design load which corresponds to the minimum excitation e(t). In this case, the mean 
load to alternating load ratio F = F m / F^j is very large, say F * 30. Time-invariant 

stiffness model prediction differs considerably from the experiment in Figure 4.16. 
Predicted amplitudes are considerably lower than the measurements, and a significant 
jump discontinuity found experimentally around the second natural frequency <Djj is not 
even predicted by this formulation. But when the periodically-varying mesh stiffness is 
considered, our model predicts the frequency response accurately including the jump 
discontinuity. The reason for a large jump around cDq for a very heavily loaded system 
(F « 30) is now clear. The second natural frequency (Ojj of the corresponding LTI 
system is nearly twice the first natural frequency (Oj. This forces the second primary 
resonance at » (Ojj to coincide with the parametric resonance at Hj, * 2coj ; 

consequently a very large jump discontinuity is developed. 

At 3/4 of the design load with F = 10, the effect of k b (t) is observed in Figure 
4. 17. Although the time-invariant model predicts the jump at (Oy, predicted amplitudes 
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Figure 4. 16 Comparison of theory with Munro’s [27] experiment at the design load. 
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Figure 4.17 Comparison of theory with Munro's [27] experiment at 3/4 of the 
design load. 
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Figure 4.18 Comparison of theory with Munro’s [27] experiment at 1/2 of 
design load. 
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Figure 4.19 Comparison of theory with Munro's [27] experiment at 1/4 of the 
design load. 
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are far below the measurements. The inclusion of kjj(t) improves predictions 
drastically. Similarly at 1/2 and 1/4 design loads, mesh stiffiiess variation as shown in 
Figures 4.18 and 4.19 respectively, affects the frequency response significantly and 
yields predictions closer to the measurements. 

4.7. CONCLUDING REMARKS 

This chapter on the non-linear dynamics of a geared rotor-bearing system with 
time-varying mesh stiffiiess k h (t), as excited by the static transmission error under a 
mean load, has resolved a number of fundamental issues. First, the interaction between 
time-varying mesh stiffness k h (t) and mean torque load has been understood. 
Second, frequency response of the corresponding LTV system has been studied, and 
the resonances associated with parametric and forced excitations have been identified. 
Third, dynamic interactions between k h (t) and system non-linearities associated with 
gear backlash and radial clearances in rolling element bearings have been investigated; a 
strong interaction between k h (t) and gear backlash is found where as the coupling 
between k^(t) and bearing non-linearities is weak. Finally, our time-varying non- 
linear formulation has yielded good predictions as compared with benchmark experiments. 
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CHAPTER V 
CONCLUSION 


5.1. SUMMARY 

In this study, linear and non-linear mathematical models of a generic geared 
rotor-bearing system shown in Figure 1.1a are developed, and several modeling issues 
which have never been addressed previously in the literature are investigated in depth. 
In Chapter I, a dynamic finite element model of the system is developed. Effects of 
several system parameters such as torsional and transverse flexibilities of the shafts and 
prime mover/load inertias are investigated, and modes of interest are identified. Three 
reduced order linear time-invariant models are developed and the conditions under 
which such models are suitable are determined by comparing the eigen-solutions with 
the finite element model results. 

In Chapter n, non-linear frequency response characteristics of a spur gear pair 
with backlash and time-invariant mesh stiffness are examined for both external and 
internal excitations. The internal excitation is of importance from the high frequency 
noise and vibration control view point and it represents the overall kinematic or static 
transmission error e(t). Such problems may be significantly different from the rattle 
problems associated with external, low frequency torque excitation. Two solution 
methods, namely the digital simulation technique and the method of harmonic balance 
have been used to develop the steady state solutions for the internal sinusoidal 
excitation. Difficulties associated with the determination of the multiple solutions at a 
given frequency in the digital simulation technique have been resolved as one must 
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search the entire initial conditions map. Such solutions and the transition frequencies 
for various impact situations are found analytically by the method of harmonic balance. 
Further, the principle of superposition can be employed to analyze the periodic 
transmission error excitation and/or combined excitation problems provided the 
excitation frequencies am sufficiently far apart from each other. Predictions compare 
reasonably well with the experimental data available in the literature. 

In Chapter ffl, non-linear frequency response characteristics of a geared rotor- 
bearing system am examined. A three degree of freedom dynamic model is developed 
which includes non-linearities associated with radial clearances in the radial rolling 
element bearings and backlash between a spur gear pair, linear time-invariant gear 
meshing stiffness is assumed. Bearing non-linear stiffness function is approximated 
for convenience sake by a simple model which is identical to that used for the gear 
mesh. This approximate bearing model has been verified by comparing the steady state 
frequency spectra. Applicability of both analytical and numerical solution techniques to 
the multi-degree of freedom non-linear problem is investigated. Proposed theory is 
validated by comparing the results with available experimental data. Several key issues 
such as non-linear modal interactions and differences between internal static 
transmission error excitation and external torque excitation are discussed. Additionally, 
parametric studies are performed to understand the effect of system parameters such as 
bearing stiffness to gear mesh stiffness ratio, alternating to mean force ratio and radial 
bearing preload to mean force ratio on the non-linear dynamic behavior. A criterion 
used to classify the steady state solutions is presented, and the conditions for chaotic, 
quasi-periodic and subharmonic steady stale solutions are determined. Two typical 
routes to chaos observed in this geared system are also identified. 
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Non-linear frequency response characteristics of a geared rotor-bearing system 
with time-varying mesh stiffness are examined in Chapter IV. The spur gear pair 
model of Chapters II and geared rotor-bearing system model of Chapter ID are 
modified to include periodic mesh stif&iess k h (t). Governing non-linear time-varying 
equations which include clearance non-linearities associated with gear backlash and 
rolling element bearings, as excited by the static transmission error e(t) under a mean 
torque load, are solved using digital simulation technique. Resonances of the 
corresponding linear time-varying (LTV) system associated with parametric and 
external excitations are identified using the method of multiple scales and digital 
simulation. Interactions between mesh stiffness variation and clearance non-linearities 
have been mvestigated; a strong interaction between time varying mesh stiffness k h (t) 
and gear backlash is found whereas the coupling between k h (t) and bearing non- 
linearities is weak. The predictions yielded by the proposed time-varying non-linear 
model agree well with the experimental results available in the literature. 

5.2 FUTURE RESEARCH AREAS 

The following topics are identified as areas of future research, based on the 
present study on non-linear dynamic analysis of geared rotor-bearing systems: 

1 . Extension of current models to study dynamics of helical and bevel gear drives, 
and multi-gear mesh systems such as planetary gear trains. 

2. Iteration procedures for the transmission error input and inclusion of the 
interactions between external and internal kinematic error excitations. Also 
inclusion of side-bands in non-linear dynamic analysis. 
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3 Design of a passive vibration control strategy using squeeze film dampers, 
investigation of semi-active and active vibration and noise control schemes 
piezoelectric actuators or hydraulic mounts. 

4 Statistical energy analysis of non-linear geared systems. 


And, 

using 
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APPENDIX A 


USER'S GUIDE FOR THE GEARED ROTOR DYNAMICS 

PROGRAM - GRD 


A.l. DESCRIPTION 

Geared Rotor Dynamics Program "GRD" is a general purpose finite element 
computer program to analyze dynamics of a system consisting of two shafts supported 
on bearings and coupled by a gear mesh. It computes natural frequencies, mode shapes 
and the vibration response of the system to static transmission error excitation, mass 
unbalances and geometric eccentricities of gears. Many rigid disks, bearings, and 
hollow shafts can be considered by GRD. Since GRD uses finite elements method, the 
system should first be discretized to small rotor elements. It is also necessary to enter 
data for each element in sequence. 

A. 2. USING THE PROGRAM 

There are two ways to input data in GRD. Data can be entered through a file or 
interactively. The first set of questions will give the user this option (the Primary 
Option) . Display on screen: 

ENTER DATA INTERACTIVELY OR FROM A DATA FILE? 

OPTION (1) : ENTER DATA INTERACTIVELY 

OPTION (2) : ENTER DATA FROM FILE 

CHOOSE OPTION 1 OR 2 (NO DEFAULT) » _ 

A.2.1. OPTION I: Interactive Data Entry 

By selecting this option, the user will be prompted to enter all the data that will be 
necessary to run this program. In the process, a new data file will be created from the 
user's input. The data entered will be reorganized in a format so that the program may 
read from it. For this reason, a new data file will be created and the user will be 
prompted to name this new file as follows: 

DATA WILL BE INPUT INTERACTIVELY, AND IN THE PROCESS, YOUR INPUT 
WILL BE ORGANIZED IN A NEW DATA FILE FOR THIS PROGRAM TO 
READ FROM. 

ENTER AN INPUT FILE NAME » 
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If an error should arise from opening this new file, the user will be given the option to 
try again: 

♦ERROR IN OPENING FILE* DO YOU WISH TO TRY AGAIN? 

TYPE (1) FOR YES OR (2) FOR NO » 

If YES, then the user will be asked the file name again. If NO, the program will 
return to the Primary Option. 


2.1.1. Data Entry Groupings 

In this section, the user will enter the necessary data. The questions will be 
grouped into the following categories. 


1. GENERAL DATA 

2. MATERIAL PROPERTIES 

3. ELEMENT PROPERTIES 

4. GEAR MESH PROPERTIES 

5. FORCED RESPONSE DATA (If user chooses this option) 

1 . General Data: 

In this section, the following general information about the overall system will be 
entered: 


1) THE NUMBER OF ROTOR ELEMENTS IN THE FIRST SHAFT 

2) THE NUMBER OF ROTOR ELEMENTS IN THE SECOND SHAFT 

3) THE NUMBER OF ROTOR ELEMENTS BEFORE THE GEAR IN THE FIRST SHAFT 

4) THE NUMBER OF ROTOR ELEMENTS BEFORE THE GEAR IN THE SECOND 
SHAFT 


As an example, consider the system shown in Figure Al. The first step is to the 
divide rotors into small pieces (finite rotor elements) to obtain a finite element model of 
the system (as illustrated in Figure A2). The shorter die rotor elements, the higher the 
accuracy and the longer the computation time. In Figure A2, each shaft is divided into 
4 pieces. The numbering can start at any end of the driving rotor and finish at any end 
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Figure Al. A Typical Geared Rotor System 
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Figure A2. Finite Element Model 
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of the driven rotor. In Figure A2 for example, we start at the left end of the first shaft 
and finished numbering at the right end of the second shaft. According to this 
configuration (Figure A2), the number of rotor elements in the first and second shafts 
are both 4 and the number of rotor elements before the pinion and gear are 2 and 6, 
respectively. 


2. Material Properties; 

In this section, material properties of the shaft considered will be entered. The 
following will be required; 


1) THE DENSITY OF THE SHAFT MATERIAL (kg/m 3 ) 

2) THE VISCOUS DAMPING COEFFICIENT (s) 

3) THE HYSTERETIC DAMPING COEFFICIENT 

4) THE MODULUS OF ELASTICITY (N/m 2 ) 

5) SHEAR MODULUS OF ELASTICITY (N/m 2 ) 


(DEFAULT =7800) 
(DEF=0) 
(DEF=0) 
(DEF=.207E12) 
(DEF=.0795E12) 


3. Element Properties; 

The main portion of data to be entered is in this section. The options given are the 
following: 


1) DEFINING THE DISK ELEMENT 

2) DEFINING THE BEARING ELEMENT 

3) DEFINING THE ROTOR ELEMENT (also the default selection) 

4) DEFINING THE GEAR ELEMENT 

5) FINISH WITH DATA INPUT FOR FIRST SHAFT 

6) FINISH WITH DATA INPUT FOR SECOND SHAFT 


The most important part in this section of questions is the order of data entry. The 
order is based on the numbered finite element model (see Figure A2). For the finite 
element model shown in Figure A2, the order of data input is as follows; 

1) DEFINING BEARING #1 

2) DEFINING ROTOR #1 

3) DEFINING ROTOR #2 

4) DEFINING GEAR#1 

5) DEFINING ROTOR #3 

6) DEFINING ROTOR #4 

7) DEFINING BEARING #2 

8) TYPE IN OPTION 5 (FINISHED WITH DATA INPUT FOR SHAFT #1 ) 

9) DEFINING BEARING #3 

10) DEFINING ROTOR #5 

1 1 ) DEFINING ROTOR #6 
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12) DEFINING GEAR#2 

13) DEFINING ROTOR #7 

14) DEFINING ROTOR #8 

15) DEFINING BEARING #4 

16) TYPE IN OPTION 6 (FINISHED WITH DATA INPUT FOR SHAFT #2) 


From the example above, one can see the ordering of data entry moving from the first 
element in shaft #1 to the last element in shaft #2. Ultimately, the program will also 
read from the data file in this order. 

According to the option selected (disk, gear, rotor element or bearing), the 
following information is required: 


The disk element - option 1: 


1 ) DISK OUTER DIAMETER IN (m) 

2) WIDTH OR THICKNESS OF DISK IN (m) 

3) MATERIAL DENSITY OF DISK IN (kg/m 3 ) (DEF=7800) 


The bearing element - option 2: 


THE STIFFNESS fN/mt 
UKxx 

2) K xy 

3) Kyx 

4) Kyy 

THE DAMPING (N-s/mt 

5) C*x 

6) C xy 

7) C yx 

8) C yy 


The rotor element - option 3: 


1 ) LENGHT OF ROTOR ELEMENT IN (m) 

2) OUTER DIAMETER OF ROTOR ELEMENT IN (m) 

3) INNER DIAMETER OF ROTOR ELEMENT IN (m) (for hollow rotor) (DEF=0) 

4) THE AXIAL LOAD IN (N) (DEF=0) 


The gear element - option 4: 


1 ) PITCH CIRCLE DIAMETER IN (m) 

2) FACE WIDTH IN (m) 
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3) GEAR MATERIAL DENSITY IN (kg/m 3 ) (DEF=7800) 

Option 5 should be chosen when the elements of the first shaft have been defined. 
The user will then be given the material properties of the second shaft., and given the 
option to make changes to this set of data. Similarly, option 6 should be entered 
when all elements of the second shaft have been entered. 


4. Gear Mesh Properties: 

The following gear mesh properties will be entered in this section: 


1) GEAR MESH STIFFNESS (N/m) 

2) GEAR MESH DAMPING (N-s/m) (DEF=0) 

3) BASE CIRCLE DIAMETER OF FIRST GEAR (m) 

4) BASE CIRCLE DIAMETER OF SECOND GEAR (m) 


Following the input of gear mesh data, the program will also ask the user for the 
number of natural frequencies wanted in the output. 


5. Forced Response Datai 

Here, the user is given the option to have forced response calculated. If the 
selection is NO, then interactive input will be complete. If the selection is YES, then 
die following options are given: 


1) WHIRLING ORBIT AT A SPECIF IED NODE 

2) DEFLECTIONS AT A SPECIFIED NODE „ _ _ _ 

3) DYNAMIC LOAD TO STATIC LOAD RATIO AT THE MESHING POINT 


Following the options menu, the general information for this section will be required: 


1) GEOMETRIC ECCENTRICITY (RUNOUT) OF GEAR 1 (m) 

2) GEOMETRIC ECCENTRICITY (RUNOUT) OF GEAR 2 (m) 

3) PEAK TO PEAK VALUE OF STATIC TRANSMISSION ERROR (m) 

4) MASS UNBALANCE OF GEAR 1 (kg-m) 

5) MASS UNBALANCE OF GEAR 2 (kg-m) 

6) NUMBER OF TEETH IN GE AR 1 (P INION) 

7) AVERAGE FORCE TRANSMITTED (N) 

8) VALUE OF MODAL DAMPING 


(DEF=0) 

(DEF-0) 

(DEF=0) 

(DEF=0) 


(DEF=.01) 


After the general information for this section has been completed, the program will 
require certain information based from the choice selected from the options menu : 
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Forced response option 1: Whirling orbit at a specified node, the following 
information will be required. 

1 ) NUMBER OF NODE AT WHICH WHIRL ORBIT IS REQUIRED 

2) ROTATIONAL SPEED OF SHAFT 1 (rad/s) 

Forced response option 2: Deflections at a specified node will be required. 

1 ) STARTING ROTATIONAL SPEED OF SHAFT 1 (rad/s) (DEF=0) 

2) UPPER LIMIT OF ROTATIONAL SPEED OF SHAFT 1 (rad/s) 

3) INCREMENT FOR ROTATIONAL SPEED (rad/s) (DEF=5) 

4) DIRECTION AT WHICH FORCED RESPONSE IS WANTED (1,2 OR 3) (DEF=1) 

5) NUMBER OF NODE AT WHICH FORCED RESPONSE IS WANTED 

In question (4), enter: 1 for deflection in pressure line direction; 2 for deflection in 
direction perpendicular to pressure line; 3 for torsional deflections 

1 ) STARTING ROTATIONAL SPEED OF SHAFT 1 (rad/s) (DEF=0) 

2) UPPER LIMIT OF ROTATIONAL SPEED OF SHAFT 1 (rad/s) 

3) INCREMENT FOR ROTATIONAL SPEED (rad/s) (DEF=5) 


Having completed this section on forced response, the program will then read from the 
new data file created from the interactive input. The output will be in an output data file 
namedFOR001.DAT . 


A.2.2. OPTION II: Enter Data From File 

With this selection, the user is then requested to give the name of the file where the 
data is stored, instead of entering data interactively. The program will then proceed to 
read from the file. If an error should occur in opening the file, then the user will be 
given the option to try again. If the user opts to try again, then the user will be asked 
the file name again. If no, then the program will return to the primary option. After the 
program has read from the file the output will be found in an output data file. 
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A. 4. SAMPLE OUTPUT FILE FOROOLDAT 


********* 






GEARED ROTOR DYNAMICS PROGRAM 
GRD 

GEAR DYNAMICS AND GEAR NOISE RESEARCH LABORATORY 
THE OHIO STATE UNIVERSITY 






THE FIRST SHAFT: 


MATERIAL PROPERTIES OF THE SHAFT: 


DENSITY OF THE MATERIAL = 

ELASTIC MODULUS = 

SHEAR MODULUS “ 

VISCOUS DAMPING COEFFICIENT.* 
HYSTERETIC LOSS FACTOR = 


0 . 78000E+04 KG/M* *3 
0 . 20700E+12 N/M**2 
0 . 79500E+11 N/M**2 
0 . OOOOOE+OO S 
O.OOOOOE+OO 


AT Z« 0.0000 M THERE EXISTS A BEARING WITH THE 
FOLLOWING STIFFNESS AND DAMPING COEFFICIENTS: 


KXX* 0.10000E+10 N/M 

KYX= O.OOOOOE+OO N/M 

CXX- 0.00000E+00 N-S/M 

CYX* 0.00000E+00 N-S/M 


KXY* O.OOOOOE+OO N/M 

KYY- 0.10000E+10 N/M 

CXY- O.OOOOOE+OO N-S/M 

CYY* 0.00000E+00 N-S/M 


AT Z- 0.0000 M THERE EXISTS A 
FOLLOWING SPECIFICATIONS: 


FINITE ROTOR ELEMENT WITH THE 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0 . 30000E-01 M 
O.OOOOOE+OO M 
0.00000E+00 N 


AT Z- 0.0500 M THERE EXISTS A FINITE ROTOR ELEMENT WITH THE 
FOLLOWING SPECIFICATIONS: 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0.50000E-01 M 
0 . 30000E-01 M 
0.00000E+00 M 
0.00000E+00 N 


AT Z- 0.1000 M THERE EXISTS A RIGID 
SPECIFICATIONS : 


DISK WITH THE FOLLOWING 


OUTER DIAMETER * 0.13500E+00 M 
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WIDTH 

MATERIAL DENSITY. . . . 


0 . 25400E-01 M 
0 . 78000E+04 KG/M**3 


AT Z— 0.1000 M THERE EXISTS A FINITE ROTOR ELEMENT WITH THE 
FOLLOWING SPECIFICATIONS: 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0 . 30000E-01 M 
0.00000E+00 M 
0 . 00000E+00 N 


AT Z- 0 
FOLLOWING 


.1500 M THERE EXISTS A FINITE ROTOR ELEMENT WITH THE 
SPECIFICATIONS: 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0.30000E-01 M 
0 . OOOOOE+OO M 
O.OOOOOE+OO N 


AT Z« 0.2000 M THERE EXISTS A BEARING WITH THE 
FOLLOWING STIFFNESS AND DAMPING COEFFICIENTS: 


KXX- 

0 . 10000E+10 

N/M 

KYX— 

0.00000E+00 

N/M 

CXX- 

0 . 00000E+00 

N-S/M 

CYX— 

0.00000E+00 

N-S/M 


THE SECOND SHAFT: 


KXY- O.OOOOOE+OO N/M 

KYY- 0 . 10000E+10 N/M 

CXY- 0.00000E+00 N-S/M 

CYY- O.OOOOOE+OO N-S/M 


MATERIAL PROPERTIES OF THE SHAFT: 


DENSITY OF THE MATERIAL - 

ELASTIC MODULUS - 

SHEAR MODULUS - 

VISCOUS DAMPING COEFFICIENT . - 
HYSTERETIC LOSS FACTOR - 


0.78000E+04 KG/M* *3 
0 . 20700E+12 N/M* *2 
0 . 79500E+11 N/M**2 
O.OOOOOE+OO S 
0.00000E+00 


AT Z- 0.0000 M THERE EXISTS A BEARING WITH THE 
FOLLOWING STIFFNESS AND DAMPING COEFFICIENTS: 


KXX- 0.10000E+10 N/M 

KYX- 0.00000E+00 N/M 

CXX- 0.00000E+00 N-S/M 

CYX- O.OOOOOE+OO N-S/M 

AT Z- 0.0000 M THERE EXISTS 
FOLLOWING SPECIFICATIONS: 


KXY- O.OOOOOE+OO N/M 

KYY- 0.10000E+10 N/M 

CXY- O.OOOOOE+OO N-S/M 

CYY- 0.00000E+00 N-S/M 

FINITE ROTOR ELEMENT WITH THE 


LENGTH OF THE ELEMENT 

. . . .« 

0 . 50000E-01 

M 

OUTER 

DIAMETER 


0. 40000E-01 
0.00000E+00 
0.00000E+00 

M 

INNER 

DIAMETER 


M 

AXIAL 

LOAD 


1*1 

N 
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AT Z« 0.0500 M THERE EXISTS A FINITE 
FOLLOWING SPECIFICATIONS: 


ROTOR ELEMENT WITH THE 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0 . 40000E-01 M 
0.00000E+00 M 
0.00000E+00 N 


AT Z- 0.1000 M THERE EXISTS A RIGID DISK WITH THE 
SPECIFICATIONS: 


FOLLOWING 


OUTER DIAMETER “ 

WIDTH - 

MATERIAL DENSITY.... - 


0 . 24500E+00 M 
0 . 25400E-01 M 
0 . 78000E+04 KG/M* *3 


AT Z« 0.1000 M THERE EXISTS A FINITE ROTOR ELEMENT 
FOLLOWING SPECIFICATIONS: 


WITH THE 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0 . 40000E-01 M 
0 . 00000E+00 M 
0 . 00000E+00 N 


AT Z- 0.1500 M THERE EXISTS A FINITE 
FOLLOWING SPECIFICATIONS: 


ROTOR ELEMENT WITH THE 


LENGTH OF THE ELEMENT 

OUTER DIAMETER 

INNER DIAMETER 

AXIAL LOAD 


0 . 50000E-01 M 
0 . 40000E-01 M 
O.OOOOOE+OO M 
0.00000E+00 N 


AT Z“ 0.2000 M THERE EXISTS A BEARING WITH THE 
FOLLOWING STIFFNESS AND DAMPING COEFFICIENTS: 


KXX* 

0 . 10000E+10 

N/M 

KYX* 

0.00000E+00 

N/M 

cxx- 

0.00000E+00 

N-S/M 

CYX- 

0.00000E+00 

N-S/M 


KXY- 

O.OOOOOE+OO 

N/M 

KYY« 

0.10000E+10 

N/M 

CXY« 

O.OOOOOE+OO 

N-S/M 

CYY- 

0.00000E+00 

N-S/M 


GEAR MESH PROPERTIES: 


BASE CIRCLE DIAMETER OF GEAR 1...- 0.1270 M 

BASE CIRCLE DIAMETER OF GEAR 2...- 0.2310 M 

AVERAGE MESH STIFFNESS “ 0.20000E+09 N/M 

AVERAGE MESH DAMPING - 0.00000E+00 N-S/M 






A. FREE VIBRATION ANALYSIS: 


FOLLOWING NATURAL FREQUENCIES AND MODESHAPES ARE CALCULATED 
FOR THE ABOVE SPECIFIED SYSTEM: 
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MODE 1 


NATURAL FREQUENCY- 0.0002 (HZ.) 

CORRESPONDING MODESHAPES : 

DISP. IN Y DISP. IN X ROT. ABOUT Y ROT. ABOUT X ROT. ABOUT Z 


- . 260E-08 0 . 000E+00 0.000E+00 -.105E-06 0.599E+01 
- . 358E-08 0 . 000E+00 0.000E+00 0.108E-06 0.599E+01 
0 . 119E-08 0.000E+00 O.OOOE+OO -.411E-07 0.599E+01 
- . 524E-08 0.000E+00 0.000E+00 -.983E-07 0.599E+01 
- . 299E-08 0.000E+00 O.OOOE+OO 0.136E-06 0.599E+01 
0 . 771E-08 0.000E+00 O.OOOE+OO 0.499E-06 0.329E+01 
0 . 253E-07 O.OOOE+OO O.OOOE+OO 0.535E-07 0.329E+01 
- . 188E-09 O.OOOE+OO O.OOOE+OO -.118E-05 0.329E+01 
- . 256E-07 O.OOOE+OO O.OOOE+OO 0.557E-07 0.329E+01 
- . 781E-08 O.OOOE+OO O.OOOE+OO 0.505E-06 0.329E+01 


MODE 2 


NATURAL FREQUENCY- 486.0262 (HZ.) 
CORRESPONDING MODESHAPES : 


IP. IN Y 

DISP. IN X 

ROT. ABOUT Y 

ROT. ABOUT X 

ROT. ABOUT Z 

0. 929E-02 

0 . 000E+00 

O.OOOE+OO 

0. 551E+01 

- . 646E+01 

0 . 261E+00 

O.OOOE+OO 

O.OOOE+OO 

0 . 410E+01 

- . 645E+0I 

0 . 375E+00 

0 . 000E+00 

O.OOOE+OO 

- . 11 9E-06 

- . 643E+01 

0 . 261E+00 

O.OOOE+OO 

O.OOOE+OO 

- . 410E+01 

- . 645E+01 

0.929E-02 

O.OOOE+OO 

O.OOOE+OO 

- . 551E+01 

- . 646E+01 

-.950E-02 

O.OOOE+OO 

O.OOOE+OO 

- . 180E+01 

0.110E+01 

918E-01 

O.OOOE+OO 

O.OOOE+OO 

-. 134E+01 

0 . 110E+01 

- . 129E+00 

O.OOOE+OO 

O.OOOE+OO 

0. 995E-09 

0 . 109E+01 

918E-01 

O.OOOE+OO 

O.OOOE+OO 

0 . 134E+01 

0.110E+01 

950E-02 

O.OOOE+OO 

O.OOOE+OO 

0.180E+01 

0 . 110E+01 
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MODE 10 


NATURAL FREQUENCY" 2479.9504 (HZ.) 

CORRESPONDING MODESHAPES: 

DISP. IN Y DISP. IN X ROT. ABOUT Y ROT. ABOUT X ROT. ABOUT Z 


0 . 124E-01 
0 . 212E+00 
0 . 279E+00 
0 . 212E+00 
0 . 124E-01 
-. 990E-02 
- . 691E-01 
- . 917E-01 
- . 691E-01 
- . 990E-02 


0 . 000E+00 
0 . 000E+00 
0.000E+00 
0.000E+00 
0.000E+00 
0 . 000E+00 
0 . 000E+00 
0 .OOOE+OO 
0. OOOE+OO 
0. OOOE+OO 


0. OOOE+OO 
0. 000E+00 
0. 000E+00 
0. 000E+00 
0. OOOE+OO 
0. OOOE+OO 
0. OOOE+OO 
0. OOOE+OO 
0 . OOOE+OO 
0 . OOOE+OO 


0 . 459E+01 
0 . 288E+01 
-.897E-10 
- . 288E+01 
- . 459E+01 
- . 134E+01 
- . 897E+00 
- . 562E-10 
0 . 897E+00 
0 . 134E+01 


0 . 978E+01 
0 . 950E+01 
0 .865E+01 
0 . 950E+01 
0.978E+01 
-. 166E+01 
- . 162E+01 
- . 147E+01 
- . 162E+01 
- . 166E+01 






B. FORCED VIBRATION ANALYSIS: 


FORCED RESPONSE AT THE SPECIFIED FREQUENCY RANGE: 

0 . 591142E-04 
0 . 675386E-04 
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APPENDIX B 


DIMENSIONAL AND NON-DIMENSIONAL QUANTITIES 
FOR THE EXPERIMENTAL TEST RIGS 


Table Bl. Dimensionless and dimensional parameters for Kubo's test rig [28]. 


DIMENSIONAL QUANTITIES 


Number of teeth 

25/25 


Igl.Ig2 [kg-m 2 (lb-ft 2 )] 

0.00115 (0.0278) 


m cl Ikg] 

0.23 


dgi,d g2 [m (in)] 

0.094 (3.7) 



rigid 


k h [N/m (lbf/in)] 

3.8E8 (2.17E6) 


F m [N (lbf)] 

2295 (1030) 


e [m (in)] 

1.92E-6 (7.56E-5) 


bh [m (in)] 

0.1 E-3 (0.0039) 



DIMENSIONLESS QUANTITIES 


C (page 44) 

0.1 

F m (page 44) 

0.06 

Fah (page 44) 

0.0192 

Fm/Fah 

3.12 
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Table B2. Dimensionless and dimensional parameters for Munro's test rig [27], 


DIMENSIONAL QUAN TITIES 

number of teeth 


32/32 



Igl,Ig2 (kg-m2 (lb-ft2)] 


0.136 (3.29) 


m gl ,m g2 [kg (lb)] 


31.1 (68.5) 


d gl’ d g2 [«*»(“)] 


0.2 (8) 



bh [m (in)] 


0.12E-3 (4.7E-3) 



Design Load (DL) 

3/4 DL 

1/2 DL 

1/4 DL 

F m [N (lbf/in FW)] 

3782 (1700) 

2836(1275) 1891 (850) 

947 (425) 

k h [N/m] 

3.44E8 

3.22E8 

3.01E8 

2.72E8 

e [m] 

3.5E-7 

1.06E-6 

1.78E-6 

2.36E-6 

DIMENSIONLESS QUANTITIES 


Design Load (DL) 

3/4 DL 

1/2 DL 

1/4 DL 

F m (equation 3.3i) 

0.183 

0.146 

0.105 

0.058 

Fata (equation 3.3j) 

0.0058 

0.0178 

0.0296 

0.0393 

K l l * k 22 (equation 3 .3f) 

0.950 

0.966 

0.983 

1.007 

k 13* k 23 (equation 3.3g) 

0.242 

0.242 

0.242 

0.242 

Cll» C22 (equation 3.3c) 

0.01 

0.01 

0.01 

0.01 

Cl3» C23 (equation 3.3d) 

0.00375 

0.00375 

0.00375 

0.00375 

£33 (equation 3.3e) 

0.015 

0.015 

0.015 

0.015 
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APPENDIX C 

LENGTH OF TRANSIENT SOLUTION AND CPU TIME 


C.l. Length of Transient Solution 

It is necessary to run the digital simulation routines for many cycles in order to 
reach the steady state solution. This can easily be detected by a time history repeating 
itself and a phase plane trajectory followed at each period of the steady state solution. 
Following are found to be the major factors which determine the length of the transient 
region of die time history: 

1. Damping Ratio: Damping ratio is inversely proportional to the time needed to 
reach steady state solution. Larger the damping, shorter the transients. Figure B1 
shows this relationship for the single degree of freedom model of gear pair. A similar 
trend is observed in Multi-degree of freedom model also. 

2. Excitation Freq uency: The transient solution is found to be longer in the 
vicinity of the resonance frequency than the off-resonance region. Figure B 1 compares 
length of transient solution versus damping ratio curves corresponding to resonance 
and off-resonance excitation frequencies. 

3. Initial Conditions: Transients are longer when the initial displacement and 
velocity are away from the mean values of the steady state solution sought. 

C.2. CPU Times 

CPU time needed in digital simulation depends on: i. the number of increments 
per period of forcing function (20-40 points per period is good enough), ii. the 
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80 



damping ratio 


Figure Cl. Length of the transient solution in a gear pair as a function of damping 
ratio for off-resonance (Q=0.5) and resonance (0=0.7) frequencies 


tolerance of the solution accuracy (lxlO' 5 to lxlO* 9 is acceptable), and iii. number of 
degrees of freedom. For the single degree of freedom model with a tolerance of lxlO 9 
and an increment of 20 points/period, nearly 10 sec. of CPU time is required to run for 
100 cycles. The three degree of freedom model with same parameters uses 1 .2 min. of 
CPU time for 100 cycles. 
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APPENDIX D 

COMPARISON OF PREDICTIONS WITH NASA MEASUREMENTS 


In this section, a comparison of the FEM predictions of the NASA Gear Noise 
Test Rig shown in Figure D1 with the measurements of NASA is given. System 
parameters and the predicted natural modes of the system have already been given in 
Tables 1.1 and 1.2, respectively. The static transmission error and mesh stiffness 
variation predictions are also shown in Figure 4.8. Here, the comparison of the 
predicted frequency response spectrum with measurements of NASA is limited to 
torsional (angular acceleration) vibrations since the set-up was not equipped for 
transverse vibration measurements. Only a qualitative agreement is observed between 
the predictions and measurements as shown in Figure D2 due to the uncertainties 
associated with the sensitivities of instruments and due to the errors involved in 
measurement such as resonance and calibration problems of angular transducers. 



gearbox 


Figure Dl. NASA Gear Noise Test Rig. 
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Figure D2. Angular acceleration spectrum of the NASA Gear Noise Test Rig at pinion 
location; a) measurements b) FEM predictions. 
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APPENDIX D 


COMPARISON OF PREDICTIONS WITH NASA MEASUREMENTS 

In this section, a comparison of the FEM predictions of the NASA Gear Noise 
Test Rig shown in Figure D1 with the measurements by NASA is given. System 
parameters and the predicted natureal modes of the system have already been given 
in Tables 1.1 and 1.2, respectively. The static transmission error and mesh stiffness 
variation predictions are also shown in Figure 4.8. Here, the comparison of the predicted 
frequency response spectrum with measurements by NASA is limited to torsional 
(angular acceleration) vibrations. Only a qualitative agreement is observed between 
the predictions and measurements as shown in Figure D2. Quantitative discrepancies 
are attributed to the resonance problems experienced with the angular transducers. 



gear box 


Figure D1 . NASA Gear Noise Test Rig. 
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